R code for discrete morphological character analyses for the paper: The patterns and modes of the evolution of disparity in Mesozoic Birds

Introduction

The R code shown here performs all of the discrete morphological character tests used in the paper “The patterns and modes of the evolution of disparity in Mesozoic Birds”. All analyses should be reproducible on your own machine as data are provided on the GitHub respository and are read in via URL rather than local paths. Please make sure you use the latest version of Claddis (installation instructions below).

Specifically, this script:

1. Employs various methods to time-scale the phylogenetic hypothes(es) used, including minimal branch length (mbl) Laurin (2004), equal Brusatte et al. (2008), and tip-dating Pyron (2011), Ronquist et al. (2012) approaches.
2. Generates pairwise distance matrices directly from the discrete characters, including using Maximum Observable Rescaled Distance (MORD; Lloyd (2016)), Generalised Euclidean Distance (GED; Wills, Briggs, and Fortey (1994)), and the Hopkins and St John (HSJ; Hopkins and St John (2018)) approaches, and]with or without inferred ancestral values.
3. Ordinates these distance matri(ces) into metric spaces using principal coordinates analysis Gower (1966).

using the 1) a distance matrix was constructed by calculating a morphological distance between each pair of taxa using the MorphDistMatrix function. We chose the and the as the distance metrics, because they are more suitable to datasets containing fossils; the distance metrics were also calculated using the ??HSJ?? approach to account for inapplicable characters. All these different approaches can be applied by changing the Arguments in MorphDistMatrix () and Claddis::ordinate_cladistic_matrix (). These different methods produced essentially same results, and thus the mbl-scaled and MORD distance methods were detailed bellow, and corresponding results were present in figures. 2) First, we explore the pattern and mode of disparity of Mesozoic birds as a whole using three metrics (sum of variances and ranges, and median distance from centroids), and then conducted intergroup comparison; 3) A nonparametric multivariate analysis of variance (PERMANOVA) test was used to test whether the three avian groups were statistically seprated in morphospace; 4) we subdivided the character matrix into six anatomical subregions (skull, vertebral column, pectoral girdle, forelimb, pelvis, and hindlimb), and conducted aforementioned comparative analyses to explore whether these subregions demonstrate different patterns; 5) finally, we performed the character exhaustion analyses.

Test reference Team (2019).

Required packages

Before we begin we need to make sure the functions used in this file are available in R by installing the required packages from CRAN.

install.packages(c("phytools", "paleotree", "devtools", "maps", "dispRity", "Claddis", "strap", "vegan", "geomorph", "devtools", "gtools", "knitr", "MASS", "phylobase", "magrittr", "kableExtra", "RColorBrewer"), dependencies = TRUE)

We also require some newer functions from the latest versions of the Claddis and dispRity packages in partciular so we will install these directly from GitHub.

devtools::install_github("graemetlloyd/Claddis")
devtools::install_github("TGuillerme/dispRity", ref = "release")

We can then load some libraries into memory, although most functions will be called with the following syntax: packagename::functionname.

library(dispRity)
library(stringr)

Loading data

Before any analyses can be perfomed the raw data must be imported into R. This includes our discrete character data, phylogenetic hypotheses and information on our taxon ages.

We will load the key raw data from the GitHub data folder into memory.

nexus.data <- Claddis::read_nexus_matrix("https://raw.githubusercontent.com/graemetlloyd/ProjectAncientBeak/main/data/matrix.nex")
tip.ages <- utils::read.table("https://raw.githubusercontent.com/graemetlloyd/ProjectAncientBeak/main/data/tipages.txt", sep = ",", header = TRUE, row.names = 1)
mpts <- ape::read.nexus("https://raw.githubusercontent.com/graemetlloyd/ProjectAncientBeak/main/data/mpts.nex")
tip.dated.tree <- ape::read.nexus("https://raw.githubusercontent.com/graemetlloyd/ProjectAncientBeak/main/data/tipdated.tre")

We can also make a strict consensus tree for our parsimony data using the consensus function in the ape package.

sc.tree <- ape::consensus(mpts)

We also need to tweak a taxon name in the tip-dated tree to make sure these all match the rest of the data.

tip.dated.tree$tip.label[tip.dated.tree$tip.label == "Mengciusornis_dentatus"] <- "Mengciusornis"

Next we can define our three main taxon groupings: Enantiornithes, Ornithuromorpha and the remaining taxa labelled simply non-Ornithothoraces). We will use these later to subset our data.

Enantiornithes <- c("Boluochia", "Concornis", "Elsornis", "Eoalulavis", "Cathayornis", "Eocathayornis", "Eoenantiornis", "Gobipteryx", "Longipteryx", "Longirostravis", "Neuquenornis", "Pengornis", "Eopengornis", "Chiappeavis", "Parapengornis", "Protopteryx", "Rapaxavis", "Shanweiniao", "Vescornis", "Piscivorenantiornis", "Linyiornis", "Sulcavis", "Bohaiornis", "Longusunguis", "Shenqiornis", "Zhouornis", "Parabohaiornis", "Fortunguavis", "Pterygornis", "Cruralispennia", "Monoenantiornis", "Shangyang", "Mirusavis", "Gretcheniao", "Dunhuangia", "Qiliania")
Ornithuromorpha <- c("Archaeorhynchus", "Schizooura", "Bellulornis", "Jianchangornis", "Songlingornis", "Longicrusavis", "Apsaravis", "Hongshanornis", "Archaeornithura", "Parahongshanornis", "Tianyuornis", "Patagopteryx", "Yixianornis", "Piscivoravis", "Iteravis", "Gansus", "Ichthyornis", "Hesperornis", "Parahesperornis", "Enaliornis", "Baptornis_advenus", "Baptornis_varneri", "Vegavis", "Mengciusornis", "Yanornis", "Similiyanornis", "Abitusavis", "Eogranivora", "Xinghaiornis", "Dingavis", "Vorona")
NonOrnithothoraces <- c("Dromaeosauridae", "Archaeopteryx", "Jeholornis", "Jinguofortis", "Chongmingia", "Sapeornis", "Confuciusornis_sanctus", "Changchengornis_hengdaoziensis", "Eoconfuciusornis_zhengi", "Confuciusornis_dui", "Yangavis")

We can now use this information to make a formal taxonGroups object for use in Claddis.

taxon_groups <- list(Enantiornithes = Enantiornithes, Ornithuromorpha = Ornithuromorpha, NonOrnithothoraces = NonOrnithothoraces)
class(taxon_groups) <- "taxonGroups"

Time-scaling phylogenetic trees

Before we can use our parsimony trees these need to be assigned branch lengths reflecting their durations and a root age “pinning” their origin in time. The tip-dated tree also needs the latter, as this is not standard output from software like BEAST, MrBayes or RevBayes (as they typically simply assume the most recent taxon will be extant, i.e., 0 Ma old).

Parsimony trees can be time-scaled in R in simple ad-hoc ways, such as setting a minimum branch length (“mbl”) of one million years.

time.trees <- lapply(mpts, function(x) paleotree::timePaleoPhy(x, tip.ages, type = "mbl", vartime = 1))
time.sc.tree <- paleotree::timePaleoPhy(sc.tree, tip.ages, type = "mbl", vartime = 1)

Alternatively, after setting every internal node as old as its’ oldest descendant taxon - a process which generates multiple problematic zero-length branches (zlbs) - these zlbs can be assigned positive duration by sharing time equally (“equal” method) with the first preceding branch of positive duration.

time.trees2 <- lapply(mpts, function(x) paleotree::timePaleoPhy(x, tip.ages, type = "equal", vartime = 1))
time.sc.tree2 <- paleotree::timePaleoPhy(sc.tree, tip.ages, type = "equal", vartime = 1)

Our tip-dated tree already has branch durations, but still needs a root age assignment. Here, as we have two extant taxa (Anas, Gallus) we can simply set this using the maximum path-length (root to Anas (or Gallus) path) in millions of years.

tip.dated.tree$root.time <- max(diag(ape::vcv(tip.dated.tree)))

Pruning crown bird taxa

As our analysis is focused on the Mesozoic, and sampling of Neornithes - the bird crown - is severely limited, we need to prune these taxa.

crown.taxa.to.remove <- c("Gallus", "Anas")
stem.time.trees <- lapply(time.trees, function(x) Claddis::drop_time_tip(x, tip_names = crown.taxa.to.remove))
stem.time.trees2 <- lapply(time.trees, function(x) Claddis::drop_time_tip(x, tip_names = crown.taxa.to.remove))
stem.time.sc.tree <- Claddis::drop_time_tip(time.sc.tree, tip_names = crown.taxa.to.remove)
stem.time.sc.tree2 <- Claddis::drop_time_tip(time.sc.tree2, tip_names = crown.taxa.to.remove)
stem.tip.dated.tree <- Claddis::drop_time_tip(tip.dated.tree, tip_names = crown.taxa.to.remove)
stem.nexus.data <- Claddis::prune_cladistic_matrix(nexus.data, taxa2prune = crown.taxa.to.remove)

Calculate discrete morphological distances and ordinate data

In order to analyse the data further morphological distances between each taxon pair need to be calculated. These distances can be used directly as a disparity measure or ordinated into a metric “morpho” space, e.g., using principal coordinates.

First we can generate a distance matrix using the Maximum Observable Rescaled Distance (“mord”) metric Lloyd (2016). This is effectively a modified version of the Gower Coefficient Gower (1971) that ensures all distances are placed on a zero to one scale. These rescale distances based on available (observable) characters as a means of dealing with missing data.)

stem.mord.dist.data <- Claddis::calculate_morphological_distances(stem.nexus.data, distance_metric = "mord", distance_transformation = "none")

Alternatively, we can use the popular Generalised Euclidean Distance (“ged”) metric. This approach deals with missing data by plugging gaps with a mean distance.

stem.ged.dist.data <- Claddis::calculate_morphological_distances(stem.nexus.data, distance_metric = "ged", distance_transformation = "none")

We can also ordinate the data into a high-dimensional space - effectively giving each taxon a set of coordinates in that space that is intended to convey its’ position relative to all other taxa in the same metric space. This enables a broader array of downstream analyses to be performed than distances alone allow.

stem.pcoa.data <- Claddis::ordinate_cladistic_matrix(stem.nexus.data, distance_metric = "mord")

The problem of inapplicable characters

An additional complication for discrete characters is the coding of some data as “inapplicable” due to hierarchical coding strategies. E.g., one character may encode the presence or absence of feathers and another the coour of those feathers. For a taxon that is coded as absent for the presence of feathers there is no logical coding for their colour. Even coding them as missing is misleading, especially if data are rescaled or an estimate for this state is made when no value is logically possible. Failing to account for this issue can lead to ranked distances that are misleading and hence distort any subsequent downstream analysis, including ordination.

A solution to this problem was put forward by Hopkins and St John (Hopkins and St John (2018); HSJ) that requires explicitly stating the hierarchical nature of character dependencies. This information is then used to encode the differences across each set of linked characters into just the character at the top of each hierarchy. We can start by stating these dependencies with a simple table.

character.dependencies <- matrix(c(35, 34, 75, 69, 173, 172), ncol = 2, byrow = TRUE, dimnames = list(c(), c("dependent_character", "independent_character")))

Next we can apply the HSJ solution using a chosen alpha value that varies the influence of the dependent characters on the distances. Here we will use 0.5 (the middle, compromise value).

stem.pcoa.data <- Claddis::ordinate_cladistic_matrix(stem.nexus.data, distance_metric = "mord", distance_inapplicable_behaviour = "hsj", character_dependencies = character.dependencies, alpha = 0.5)

Incalculable distances and removing poorly known taxa

Due to missing data some pairwise distances can be incalculable. For example, if one fosil taxon is only know from cranial remains and another from only postcranial remains there are no charaters that can be coded for both taxa and hence no way to calculate a distance, regardless of wehther mord or ged is used. If applying distance-only metrics this is not a problem, but if we desire to ordinate our data into a metric space there must be a value for every pairwise distance and hence Claddis will iteratively remove the taxon that cuses the most incalculable distances until the matrix is complete.

We can find any removed taxa we found using the appropriate output.

stem.pcoa.data$removed_taxa

For consistency we now need to remove these taxa from our phylogenetic hypotheses. Note that we do this after they have contributed information to time-scaling those trees as this is best practice for optimal branch-length estimation. In other words, we should use the maximum available information for every test until it is insufficent.

stem.time.trees <- lapply(stem.time.trees, function(x) Claddis::drop_time_tip(x, tip_names = stem.pcoa.data$removed_taxa))
stem.time.trees2 <- lapply(stem.time.trees2, function(x) Claddis::drop_time_tip(x, tip_names = stem.pcoa.data$removed_taxa))
stem.time.sc.tree <- Claddis::drop_time_tip(stem.time.sc.tree, tip_names = stem.pcoa.data$removed_taxa)
stem.time.sc.tree2 <- Claddis::drop_time_tip(stem.time.sc.tree2, tip_names = stem.pcoa.data$removed_taxa)
stem.tip.dated.tree <- Claddis::drop_time_tip(stem.tip.dated.tree, tip_names = stem.pcoa.data$removed_taxa)
Enantiornithes <- setdiff(Enantiornithes, stem.pcoa.data$removed_taxa)
Ornithuromorpha <- setdiff(Ornithuromorpha, stem.pcoa.data$removed_taxa)
NonOrnithothoraces <- setdiff(NonOrnithothoraces, stem.pcoa.data$removed_taxa)
taxon_groups <- lapply(taxon_groups, function(x) {indices.to.remove <- sort(match(stem.pcoa.data$removed_taxa, x)); if(length(indices.to.remove) > 0) x <- x[-indices.to.remove]; x})
class(taxon_groups) <- "taxonGroups"

Plotting our discrete character morphospaces

As with any analysis it can be important to visualise the data to identify general patterns and guide interpretation. Here we will use some of the plotting function froms Claddis.

First we will do some simple bivariate plots for each combination of the first three axes.

Claddis::plot_morphospace(stem.pcoa.data, x_axis = 1, y_axis = 2, plot_taxon_names = TRUE, taxon_groups = taxon_groups, plot_convex_hulls = TRUE)

Claddis::plot_morphospace(stem.pcoa.data, x_axis = 1, y_axis = 3, plot_taxon_names = TRUE, taxon_groups = taxon_groups, plot_convex_hulls = TRUE)

Claddis::plot_morphospace(stem.pcoa.data, x_axis = 2, y_axis = 3, plot_taxon_names = TRUE, taxon_groups = taxon_groups, plot_convex_hulls = TRUE)

As we can see from the axis labels an issue with discrete character data sets is that their metric spaces tend to distribute the variance over a large number of axes, making it difficult to easily visualise the full range of variance. We can thus llok at more axes with the plot_multi_morphospace function.

Claddis::plot_multi_morphospace(stem.pcoa.data, n_axes = 8, plot_taxon_names = FALSE, taxon_groups = taxon_groups, plot_convex_hulls = TRUE)

We can also make a stacked morphospace plot Foote (1993), but first we need to make a new tip ages object.

pcoa.ages.data <- tip.ages[rownames(stem.pcoa.data$distance_matrix), ]
colnames(pcoa.ages.data) <- c("fad", "lad")

As well as a timeBins object that defines the slices we want to plot.

time_bins <- matrix(c(174.0, 145.0, 145.0, 125.0, 125.0, 100.5, 100.5, 66.0), ncol = 2, byrow = TRUE, dimnames = list(c("Upper Jurassic", "Berriasian-Barremian", "Aptian-Albian", "Upper Cretaceous"), c("fad", "lad")))
class(time_bins) <- "timeBins"

Now we can make the stack plot.

Claddis::plot_morphospace_stack(stem.pcoa.data, taxon_ages = pcoa.ages.data, taxon_groups, time_bins, palette = "viridis", group_legend_position = "bottom_right")

Calculate disparity metrics

We can also use our data to calculate disparity metrics, either for taxonomic groups or time bins.

A simple ordination-free metric is simply the avergae (mean) pairwsie distance between each taxon in our three groups.

Claddis::calculate_MPD(stem.mord.dist.data, taxon_groups)

Alternatively we can weight each distance by the amount of available (observable) information and take the mean of that.

Claddis::calculate_WMPD(stem.mord.dist.data, taxon_groups)

We can also apply the same measures for our time bins instead.

Claddis::calculate_MPD(stem.mord.dist.data, taxon_groups = Claddis::assign_taxa_to_bins(taxon_ages = pcoa.ages.data, time_bins))
Claddis::calculate_WMPD(stem.mord.dist.data, taxon_groups = Claddis::assign_taxa_to_bins(taxon_ages = pcoa.ages.data, time_bins))

Many disparity metrics can also be generated from our ordinations apce, such as sum of ranges…

c(Enantiornithes = sum(abs(apply(apply(stem.pcoa.data$vectors[Enantiornithes, ], 2, range), 2, diff))), Ornithuromorpha = sum(abs(apply(apply(stem.pcoa.data$vectors[Ornithuromorpha, ], 2, range), 2, diff))), NonOrnithothoraces = sum(abs(apply(apply(stem.pcoa.data$vectors[NonOrnithothoraces, ], 2, range), 2, diff))))

…sum of variances…

c(Enantiornithes = sum(apply(stem.pcoa.data$vectors[Enantiornithes, ], 2, var)), Ornithuromorpha = sum(apply(stem.pcoa.data$vectors[Ornithuromorpha, ], 2, var)), NonOrnithothoraces = sum(apply(stem.pcoa.data$vectors[NonOrnithothoraces, ], 2, var)))

…product of ranges….

c(Enantiornithes = prod(abs(apply(apply(stem.pcoa.data$vectors[Enantiornithes, ], 2, range), 2, diff))), Ornithuromorpha = prod(abs(apply(apply(stem.pcoa.data$vectors[Ornithuromorpha, ], 2, range), 2, diff))), NonOrnithothoraces = prod(abs(apply(apply(stem.pcoa.data$vectors[NonOrnithothoraces, ], 2, range), 2, diff))))

…or product of variances.

c(Enantiornithes = prod(apply(stem.pcoa.data$vectors[Enantiornithes, ], 2, var)), Ornithuromorpha = prod(apply(stem.pcoa.data$vectors[Ornithuromorpha, ], 2, var)), NonOrnithothoraces = prod(apply(stem.pcoa.data$vectors[NonOrnithothoraces, ], 2, var)))

We can also make a simple set of barplots to compare the results of these metrics.

par(mfrow = c(3, 2))
barplot(Claddis::calculate_MPD(stem.mord.dist.data, taxon_groups), main = "Mean Pairwise Distance")
barplot(Claddis::calculate_WMPD(stem.mord.dist.data, taxon_groups), main = "Weighted Mean Pairwise Distance")
barplot(c(Enantiornithes = sum(abs(apply(apply(stem.pcoa.data$vectors[Enantiornithes, ], 2, range), 2, diff))), Ornithuromorpha = sum(abs(apply(apply(stem.pcoa.data$vectors[Ornithuromorpha, ], 2, range), 2, diff))), NonOrnithothoraces = sum(abs(apply(apply(stem.pcoa.data$vectors[NonOrnithothoraces, ], 2, range), 2, diff)))), main = "Sum Of Ranges")
barplot(c(Enantiornithes = sum(apply(stem.pcoa.data$vectors[Enantiornithes, ], 2, var)), Ornithuromorpha = sum(apply(stem.pcoa.data$vectors[Ornithuromorpha, ], 2, var)), NonOrnithothoraces = sum(apply(stem.pcoa.data$vectors[NonOrnithothoraces, ], 2, var))), main = "Sum Of Variances")
barplot(c(Enantiornithes = prod(abs(apply(apply(stem.pcoa.data$vectors[Enantiornithes, ], 2, range), 2, diff))), Ornithuromorpha = prod(abs(apply(apply(stem.pcoa.data$vectors[Ornithuromorpha, ], 2, range), 2, diff))), NonOrnithothoraces = prod(abs(apply(apply(stem.pcoa.data$vectors[NonOrnithothoraces, ], 2, range), 2, diff)))), main = "Product Of Ranges")
barplot(c(Enantiornithes = prod(apply(stem.pcoa.data$vectors[Enantiornithes, ], 2, var)), Ornithuromorpha = prod(apply(stem.pcoa.data$vectors[Ornithuromorpha, ], 2, var)), NonOrnithothoraces = prod(apply(stem.pcoa.data$vectors[NonOrnithothoraces, ], 2, var))), main = "Product Of Variances")

Note that products can be problematic when they involve multiplying through by the very low variance high axes. We can thus trim our data to just the first 30 axes.

pcoa.data <- as.data.frame(stem.pcoa.data$vectors[, 1:30])
pcoa.data <- data.matrix(pcoa.data)

And as we are going to use the dispRity package convert our ages to the preferred format.

disprity.ages.data <- pcoa.ages.data
colnames(disprity.ages.data) <- c("FAD", "LAD")

And set some new time bins for the time series analyses.

time.bins <- c(170, 145, 129.4, 100.5, 86.3, 66)

We can now use our various phylogenetc hypotheses and time bin scheme to create dispRity objects.

stem.time.sc.trees.binned.morphospace <- lapply(stem.time.trees, function(x) stem.time.sc.tree.binned.morphospace <- dispRity::chrono.subsets(data = pcoa.data, tree = x, method = "discrete", time = time.bins, inc.nodes = FALSE, FADLAD = disprity.ages.data))
stem.time.sc.trees2.binned.morphospace <- lapply(stem.time.trees2, function(x) stem.time.sc.tree2.binned.morphospace <- dispRity::chrono.subsets(data = pcoa.data, tree = x, method = "discrete", time = time.bins, inc.nodes = FALSE, FADLAD = disprity.ages.data))
stem.time.sc.tree.binned.morphospace <- dispRity::chrono.subsets(data = pcoa.data, tree = stem.time.sc.tree, method = "discrete", time = time.bins, inc.nodes = FALSE, FADLAD = disprity.ages.data)
stem.time.sc.tree2.binned.morphospace <- dispRity::chrono.subsets(data = pcoa.data, tree = stem.time.sc.tree2, method = "discrete", time = time.bins, inc.nodes = FALSE, FADLAD = disprity.ages.data)
stem.tip.dated.tree.binned.morphospace <- dispRity::chrono.subsets(data = pcoa.data, tree = stem.tip.dated.tree, method = "discrete", time = time.bins, inc.nodes = FALSE, FADLAD = disprity.ages.data)

We can also use dispRity to bootstrap our data and give a better sense of our uncertianty.

stem.time.sc.trees.boot.bin.morphospace <- lapply(stem.time.sc.trees.binned.morphospace, function(x) dispRity::boot.matrix(x, bootstrap = 1000))
stem.time.sc.trees2.boot.bin.morphospace <- lapply(stem.time.sc.trees2.binned.morphospace, function(x) dispRity::boot.matrix(x, bootstrap = 1000))
stem.time.sc.tree.boot.bin.morphospace <- dispRity::boot.matrix(stem.time.sc.tree.binned.morphospace, bootstrap = 1000)
stem.time.sc.tree2.boot.bin.morphospace <- dispRity::boot.matrix(stem.time.sc.tree2.binned.morphospace, bootstrap = 1000)
stem.tip.dated.tree.boot.bin.morphospace <- dispRity::boot.matrix(stem.tip.dated.tree.binned.morphospace, bootstrap = 1000)

Now we can calculate some disparity metrics again, such as sum of variances.

stem.time.sc.trees.boot.bin.morphospace.sov <- lapply(stem.time.sc.trees.boot.bin.morphospace, function(x) dispRity::dispRity(x, metric = c(sum, variances)))
stem.time.sc.trees2.boot.bin.morphospace.sov <- lapply(stem.time.sc.trees2.boot.bin.morphospace, function(x) dispRity::dispRity(x, metric = c(sum, variances)))
stem.time.sc.tree.boot.bin.morphospace.sov <- dispRity::dispRity(stem.time.sc.tree.boot.bin.morphospace, metric = c(sum, variances))
stem.time.sc.tree2.boot.bin.morphospace.sov <- dispRity::dispRity(stem.time.sc.tree2.boot.bin.morphospace, metric = c(sum, variances))
stem.tip.dated.tree.boot.bin.morphospace.sov <- dispRity::dispRity(stem.tip.dated.tree.boot.bin.morphospace, metric = c(sum, variances))

And visualise them.

plot(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = c(0, 0, 0, 0, 0), xlim = c(170, 66), ylim = c(0, 4.5), type = "n", xlab = "Time (Ma)", ylab = "Sum Of Variance")
trees1 <- lapply(stem.time.sc.trees.boot.bin.morphospace.sov, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.1), lwd = 0.3)))
trees2 <- lapply(stem.time.sc.trees2.boot.bin.morphospace.sov, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.1), lwd = 0.3)))
trees3 <- apply(do.call(rbind, lapply(stem.time.sc.tree.boot.bin.morphospace.sov$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.2), lwd = 0.3))
trees4 <- apply(do.call(rbind, lapply(stem.time.sc.tree2.boot.bin.morphospace.sov$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.2), lwd = 0.3))
trees5 <- apply(do.call(rbind, lapply(stem.tip.dated.tree.boot.bin.morphospace.sov$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = apply(do.call(rbind, lapply(stem.time.sc.trees.boot.bin.morphospace.sov, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = apply(do.call(rbind, lapply(stem.time.sc.trees.boot.bin.morphospace.sov, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "red", lwd = 8)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = apply(do.call(rbind, lapply(stem.time.sc.trees2.boot.bin.morphospace.sov, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = apply(do.call(rbind, lapply(stem.time.sc.trees2.boot.bin.morphospace.sov, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "blue", lwd = 8)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.time.sc.tree.boot.bin.morphospace.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.time.sc.tree.boot.bin.morphospace.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "red", lwd = 8)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.time.sc.tree2.boot.bin.morphospace.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.time.sc.tree2.boot.bin.morphospace.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "blue", lwd = 8)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.tip.dated.tree.boot.bin.morphospace.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.tip.dated.tree.boot.bin.morphospace.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)

Or sum of ranges.

stem.time.sc.trees.boot.bin.morphospace.sor <- lapply(stem.time.sc.trees.boot.bin.morphospace, function(x) dispRity::dispRity(x, metric = c(sum, ranges)))
stem.time.sc.trees2.boot.bin.morphospace.sor <- lapply(stem.time.sc.trees2.boot.bin.morphospace, function(x) dispRity::dispRity(x, metric = c(sum, ranges)))
stem.time.sc.tree.boot.bin.morphospace.sor <- dispRity::dispRity(stem.time.sc.tree.boot.bin.morphospace, metric = c(sum, ranges))
stem.time.sc.tree2.boot.bin.morphospace.sor <- dispRity::dispRity(stem.time.sc.tree2.boot.bin.morphospace, metric = c(sum, ranges))
stem.tip.dated.tree.boot.bin.morphospace.sor <- dispRity::dispRity(stem.tip.dated.tree.boot.bin.morphospace, metric = c(sum, ranges))

And visualise them.

plot(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = c(0, 0, 0, 0, 0), xlim = c(170, 66), ylim = c(0, 45), type = "n", xlab = "Time (Ma)", ylab = "Sum Of Ranges")
trees1 <- lapply(stem.time.sc.trees.boot.bin.morphospace.sor, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.1), lwd = 0.3)))
trees2 <- lapply(stem.time.sc.trees2.boot.bin.morphospace.sor, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.1), lwd = 0.3)))
trees3 <- apply(do.call(rbind, lapply(stem.time.sc.tree.boot.bin.morphospace.sor$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.2), lwd = 0.3))
trees4 <- apply(do.call(rbind, lapply(stem.time.sc.tree2.boot.bin.morphospace.sor$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.2), lwd = 0.3))
trees5 <- apply(do.call(rbind, lapply(stem.tip.dated.tree.boot.bin.morphospace.sor$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = apply(do.call(rbind, lapply(stem.time.sc.trees.boot.bin.morphospace.sor, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = apply(do.call(rbind, lapply(stem.time.sc.trees.boot.bin.morphospace.sor, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "red", lwd = 8)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = apply(do.call(rbind, lapply(stem.time.sc.trees2.boot.bin.morphospace.sor, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = apply(do.call(rbind, lapply(stem.time.sc.trees2.boot.bin.morphospace.sor, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "blue", lwd = 8)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.time.sc.tree.boot.bin.morphospace.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.time.sc.tree.boot.bin.morphospace.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "red", lwd = 8)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.time.sc.tree2.boot.bin.morphospace.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.time.sc.tree2.boot.bin.morphospace.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "blue", lwd = 8)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.tip.dated.tree.boot.bin.morphospace.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.tip.dated.tree.boot.bin.morphospace.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)

Or centroid medians.

stem.time.sc.trees.boot.bin.morphospace.moc <- lapply(stem.time.sc.trees.boot.bin.morphospace, function(x) dispRity::dispRity(x, metric = c(median, centroids)))
stem.time.sc.trees2.boot.bin.morphospace.moc <- lapply(stem.time.sc.trees2.boot.bin.morphospace, function(x) dispRity::dispRity(x, metric = c(median, centroids)))
stem.time.sc.tree.boot.bin.morphospace.moc <- dispRity::dispRity(stem.time.sc.tree.boot.bin.morphospace, metric = c(median, centroids))
stem.time.sc.tree2.boot.bin.morphospace.moc <- dispRity::dispRity(stem.time.sc.tree2.boot.bin.morphospace, metric = c(median, centroids))
stem.tip.dated.tree.boot.bin.morphospace.moc <- dispRity::dispRity(stem.tip.dated.tree.boot.bin.morphospace, metric = c(median, centroids))

And visualise them.

plot(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = c(0, 0, 0, 0, 0), xlim = c(170, 66), ylim = c(0, 2.1), type = "n", xlab = "Time (Ma)", ylab = "Centroid Median")
trees1 <- lapply(stem.time.sc.trees.boot.bin.morphospace.moc, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.1), lwd = 0.3)))
trees2 <- lapply(stem.time.sc.trees2.boot.bin.morphospace.moc, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.1), lwd = 0.3)))
trees3 <- apply(do.call(rbind, lapply(stem.time.sc.tree.boot.bin.morphospace.moc$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.2), lwd = 0.3))
trees4 <- apply(do.call(rbind, lapply(stem.time.sc.tree2.boot.bin.morphospace.moc$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.2), lwd = 0.3))
trees5 <- apply(do.call(rbind, lapply(stem.tip.dated.tree.boot.bin.morphospace.moc$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = apply(do.call(rbind, lapply(stem.time.sc.trees.boot.bin.morphospace.moc, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = apply(do.call(rbind, lapply(stem.time.sc.trees.boot.bin.morphospace.moc, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "red", lwd = 8)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = apply(do.call(rbind, lapply(stem.time.sc.trees2.boot.bin.morphospace.moc, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = apply(do.call(rbind, lapply(stem.time.sc.trees2.boot.bin.morphospace.moc, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "blue", lwd = 8)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.time.sc.tree.boot.bin.morphospace.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.time.sc.tree.boot.bin.morphospace.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "red", lwd = 8)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.time.sc.tree2.boot.bin.morphospace.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.time.sc.tree2.boot.bin.morphospace.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "blue", lwd = 8)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.tip.dated.tree.boot.bin.morphospace.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.tip.dated.tree.boot.bin.morphospace.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)

We can also subset the data by taxonomic group.

stem.aves.groups <- dispRity::custom.subsets(pcoa.data, group = list(NonOrnithothoraces = match(NonOrnithothoraces, rownames(pcoa.data)), Enantiornithes = match(Enantiornithes, rownames(pcoa.data)), Ornithuromorpha = match(Ornithuromorpha, rownames(pcoa.data))))

Bootstrap it.

stem.aves.groups.bootstrap <- dispRity::boot.matrix(stem.aves.groups, bootstrap = 1000)

And calculate the same disparity metrics as the time series.

stem.aves.groups.bootstrap.sov <- dispRity::dispRity(stem.aves.groups.bootstrap, metric = c(sum, variances))
stem.aves.groups.bootstrap.sor <- dispRity::dispRity(stem.aves.groups.bootstrap, metric = c(sum, ranges))
stem.aves.groups.bootstrap.moc <- dispRity::dispRity(stem.aves.groups.bootstrap, metric = c(median, centroids))

And visualise the results.

boxplot(lapply(stem.aves.groups.bootstrap.sov$disparity, function(x) unlist(x)), ylab = "Sum Of Variance")

boxplot(lapply(stem.aves.groups.bootstrap.sor$disparity, function(x) unlist(x)), ylab = "Sum Of Ranges")

boxplot(lapply(stem.aves.groups.bootstrap.moc$disparity, function(x) unlist(x)), ylab = "Centroid Median")

Enantiornithes versus Ornithuromorpha comparisons

Our main aim is to compare the morphospaces of the two major sister-clades of Mesozoic birds, Enantiornithes and Ornithuromorpha. We will begin by creating separate PCoA data sets.

Enantiornithes.pcoa.data <- pcoa.data[Enantiornithes, ]
Ornithuromorpha.pcoa.data <- pcoa.data[Ornithuromorpha, ]

We can also create Enantiornithes-only trees.

Enantiornithes.stem.time.trees <- lapply(stem.time.trees, function(x) Claddis::drop_time_tip(x, tip_names = setdiff(x$tip.label, Enantiornithes)))
Enantiornithes.stem.time.trees2 <- lapply(stem.time.trees2, function(x) Claddis::drop_time_tip(x, tip_names = setdiff(x$tip.label, Enantiornithes)))
Enantiornithes.stem.time.sc.tree <- Claddis::drop_time_tip(stem.time.sc.tree, tip_names = setdiff(stem.time.sc.tree$tip.label, Enantiornithes))
Enantiornithes.stem.time.sc.tree2 <- Claddis::drop_time_tip(stem.time.sc.tree2, tip_names = setdiff(stem.time.sc.tree2$tip.label, Enantiornithes))
Enantiornithes.stem.tip.dated.tree <- Claddis::drop_time_tip(stem.tip.dated.tree, tip_names = setdiff(stem.tip.dated.tree$tip.label, Enantiornithes))

And Ornithuromorpha-only trees.

Ornithuromorpha.stem.time.trees <- lapply(stem.time.trees, function(x) Claddis::drop_time_tip(x, tip_names = setdiff(x$tip.label, Ornithuromorpha)))
Ornithuromorpha.stem.time.trees2 <- lapply(stem.time.trees2, function(x) Claddis::drop_time_tip(x, tip_names = setdiff(x$tip.label, Ornithuromorpha)))
Ornithuromorpha.stem.time.sc.tree <- Claddis::drop_time_tip(stem.time.sc.tree, tip_names = setdiff(stem.time.sc.tree$tip.label, Ornithuromorpha))
Ornithuromorpha.stem.time.sc.tree2 <- Claddis::drop_time_tip(stem.time.sc.tree2, tip_names = setdiff(stem.time.sc.tree2$tip.label, Ornithuromorpha))
Ornithuromorpha.stem.tip.dated.tree <- Claddis::drop_time_tip(stem.tip.dated.tree, tip_names = setdiff(stem.tip.dated.tree$tip.label, Ornithuromorpha))

We can also create some new time bins that only span our two major clades.

comparison_time_bins <- c(150, 125, 100.5, 66)

Now we can build new dispRity objects for Enantiornithes.

Enantiornithes.trees.binned.morphospace <- lapply(Enantiornithes.stem.time.trees, function(x) dispRity::chrono.subsets(data = Enantiornithes.pcoa.data, tree = x, method = "discrete", time = comparison_time_bins, inc.nodes = FALSE, FADLAD = disprity.ages.data[Enantiornithes, ]))
Enantiornithes.trees2.binned.morphospace <- lapply(Enantiornithes.stem.time.trees2, function(x) dispRity::chrono.subsets(data = Enantiornithes.pcoa.data, tree = x, method = "discrete", time = comparison_time_bins, inc.nodes = FALSE, FADLAD = disprity.ages.data[Enantiornithes, ]))
Enantiornithes.sc.tree.binned.morphospace <- dispRity::chrono.subsets(data = Enantiornithes.pcoa.data, tree = Enantiornithes.stem.time.sc.tree, method = "discrete", time = comparison_time_bins, inc.nodes = FALSE, FADLAD = disprity.ages.data[Enantiornithes, ])
Enantiornithes.sc.tree2.binned.morphospace <- dispRity::chrono.subsets(data = Enantiornithes.pcoa.data, tree = Enantiornithes.stem.time.sc.tree2, method = "discrete", time = comparison_time_bins, inc.nodes = FALSE, FADLAD = disprity.ages.data[Enantiornithes, ])
Enantiornithes.tip.dated.tree.binned.morphospace <- dispRity::chrono.subsets(data = Enantiornithes.pcoa.data, tree = Enantiornithes.stem.tip.dated.tree, method = "discrete", time = comparison_time_bins, inc.nodes = FALSE, FADLAD = disprity.ages.data[Enantiornithes, ])

And Ornithuromorpha.

Ornithuromorpha.trees.binned.morphospace <- lapply(Ornithuromorpha.stem.time.trees, function(x) dispRity::chrono.subsets(data = Ornithuromorpha.pcoa.data, tree = x, method = "discrete", time = comparison_time_bins, inc.nodes = FALSE, FADLAD = disprity.ages.data[Ornithuromorpha, ]))
Ornithuromorpha.trees2.binned.morphospace <- lapply(Ornithuromorpha.stem.time.trees2, function(x) dispRity::chrono.subsets(data = Ornithuromorpha.pcoa.data, tree = x, method = "discrete", time = comparison_time_bins, inc.nodes = FALSE, FADLAD = disprity.ages.data[Ornithuromorpha, ]))
Ornithuromorpha.sc.tree.binned.morphospace <- dispRity::chrono.subsets(data = Ornithuromorpha.pcoa.data, tree = Ornithuromorpha.stem.time.sc.tree, method = "discrete", time = comparison_time_bins, inc.nodes = FALSE, FADLAD = disprity.ages.data[Ornithuromorpha, ])
Ornithuromorpha.sc.tree2.binned.morphospace <- dispRity::chrono.subsets(data = Ornithuromorpha.pcoa.data, tree = Ornithuromorpha.stem.time.sc.tree2, method = "discrete", time = comparison_time_bins, inc.nodes = FALSE, FADLAD = disprity.ages.data[Ornithuromorpha, ])
Ornithuromorpha.tip.dated.tree.binned.morphospace <- dispRity::chrono.subsets(data = Ornithuromorpha.pcoa.data, tree = Ornithuromorpha.stem.tip.dated.tree, method = "discrete", time = comparison_time_bins, inc.nodes = FALSE, FADLAD = disprity.ages.data[Ornithuromorpha, ])

And bootstrap the data.

Enantiornithes.boot.trees.binned.morphospace <- lapply(Enantiornithes.trees.binned.morphospace, function(x) dispRity::boot.matrix(x, bootstrap = 1000))
Enantiornithes.boot.trees2.binned.morphospace <- lapply(Enantiornithes.trees2.binned.morphospace, function(x) dispRity::boot.matrix(x, bootstrap = 1000))
Enantiornithes.boot.sc.tree.binned.morphospace <- dispRity::boot.matrix(Enantiornithes.sc.tree.binned.morphospace, bootstrap = 1000)
Enantiornithes.boot.sc.tree2.binned.morphospace <- dispRity::boot.matrix(Enantiornithes.sc.tree2.binned.morphospace, bootstrap = 1000)
Enantiornithes.boot.tip.dated.tree.binned.morphospace <- dispRity::boot.matrix(Enantiornithes.tip.dated.tree.binned.morphospace, bootstrap = 1000)
Ornithuromorpha.boot.trees.binned.morphospace <- lapply(Ornithuromorpha.trees.binned.morphospace, function(x) dispRity::boot.matrix(x, bootstrap = 1000))
Ornithuromorpha.boot.trees2.binned.morphospace <- lapply(Ornithuromorpha.trees2.binned.morphospace, function(x) dispRity::boot.matrix(x, bootstrap = 1000))
Ornithuromorpha.boot.sc.tree.binned.morphospace <- dispRity::boot.matrix(Ornithuromorpha.sc.tree.binned.morphospace, bootstrap = 1000)
Ornithuromorpha.boot.sc.tree2.binned.morphospace <- dispRity::boot.matrix(Ornithuromorpha.sc.tree2.binned.morphospace, bootstrap = 1000)
Ornithuromorpha.boot.tip.dated.tree.binned.morphospace <- dispRity::boot.matrix(Ornithuromorpha.tip.dated.tree.binned.morphospace, bootstrap = 1000)

Calculate sum of variances.

Enantiornithes.boot.trees.disparity.sov <- lapply(Enantiornithes.boot.trees.binned.morphospace, function(x) dispRity::dispRity(x, metric = c(sum, variances)))
Enantiornithes.boot.trees2.disparity.sov <- lapply(Enantiornithes.boot.trees2.binned.morphospace, function(x) dispRity::dispRity(x, metric = c(sum, variances)))
Enantiornithes.boot.sc.tree.disparity.sov <- dispRity::dispRity(Enantiornithes.boot.sc.tree.binned.morphospace, metric = c(sum, variances))
Enantiornithes.boot.sc.tree2.disparity.sov <- dispRity::dispRity(Enantiornithes.boot.sc.tree2.binned.morphospace, metric = c(sum, variances))
Enantiornithes.boot.tip.dated.tree.disparity.sov <- dispRity::dispRity(Enantiornithes.boot.tip.dated.tree.binned.morphospace, metric = c(sum, variances))
Ornithuromorpha.boot.trees.disparity.sov <- lapply(Ornithuromorpha.boot.trees.binned.morphospace, function(x) dispRity::dispRity(x, metric = c(sum, variances)))
Ornithuromorpha.boot.trees2.disparity.sov <- lapply(Ornithuromorpha.boot.trees2.binned.morphospace, function(x) dispRity::dispRity(x, metric = c(sum, variances)))
Ornithuromorpha.boot.sc.tree.disparity.sov <- dispRity::dispRity(Ornithuromorpha.boot.sc.tree.binned.morphospace, metric = c(sum, variances))
Ornithuromorpha.boot.sc.tree2.disparity.sov <- dispRity::dispRity(Ornithuromorpha.boot.sc.tree2.binned.morphospace, metric = c(sum, variances))
Ornithuromorpha.boot.tip.dated.tree.disparity.sov <- dispRity::dispRity(Ornithuromorpha.boot.tip.dated.tree.binned.morphospace, metric = c(sum, variances))

Calculate sum of ranges.

Enantiornithes.boot.trees.disparity.sor <- lapply(Enantiornithes.boot.trees.binned.morphospace, function(x) dispRity::dispRity(x, metric = c(sum, ranges)))
Enantiornithes.boot.trees2.disparity.sor <- lapply(Enantiornithes.boot.trees2.binned.morphospace, function(x) dispRity::dispRity(x, metric = c(sum, ranges)))
Enantiornithes.boot.sc.tree.disparity.sor <- dispRity::dispRity(Enantiornithes.boot.sc.tree.binned.morphospace, metric = c(sum, ranges))
Enantiornithes.boot.sc.tree2.disparity.sor <- dispRity::dispRity(Enantiornithes.boot.sc.tree2.binned.morphospace, metric = c(sum, ranges))
Enantiornithes.boot.tip.dated.tree.disparity.sor <- dispRity::dispRity(Enantiornithes.boot.tip.dated.tree.binned.morphospace, metric = c(sum, ranges))
Ornithuromorpha.boot.trees.disparity.sor <- lapply(Ornithuromorpha.boot.trees.binned.morphospace, function(x) dispRity::dispRity(x, metric = c(sum, ranges)))
Ornithuromorpha.boot.trees2.disparity.sor <- lapply(Ornithuromorpha.boot.trees2.binned.morphospace, function(x) dispRity::dispRity(x, metric = c(sum, ranges)))
Ornithuromorpha.boot.sc.tree.disparity.sor <- dispRity::dispRity(Ornithuromorpha.boot.sc.tree.binned.morphospace, metric = c(sum, ranges))
Ornithuromorpha.boot.sc.tree2.disparity.sor <- dispRity::dispRity(Ornithuromorpha.boot.sc.tree2.binned.morphospace, metric = c(sum, ranges))
Ornithuromorpha.boot.tip.dated.tree.disparity.sor <- dispRity::dispRity(Ornithuromorpha.boot.tip.dated.tree.binned.morphospace, metric = c(sum, ranges))

Calculate median centroids.

Enantiornithes.boot.trees.disparity.moc <- lapply(Enantiornithes.boot.trees.binned.morphospace, function(x) dispRity::dispRity(x, metric = c(median, centroids)))
Enantiornithes.boot.trees2.disparity.moc <- lapply(Enantiornithes.boot.trees2.binned.morphospace, function(x) dispRity::dispRity(x, metric = c(median, centroids)))
Enantiornithes.boot.sc.tree.disparity.moc <- dispRity::dispRity(Enantiornithes.boot.sc.tree.binned.morphospace, metric = c(median, centroids))
Enantiornithes.boot.sc.tree2.disparity.moc <- dispRity::dispRity(Enantiornithes.boot.sc.tree2.binned.morphospace, metric = c(median, centroids))
Enantiornithes.boot.tip.dated.tree.disparity.moc <- dispRity::dispRity(Enantiornithes.boot.tip.dated.tree.binned.morphospace, metric = c(median, centroids))
Ornithuromorpha.boot.trees.disparity.moc <- lapply(Ornithuromorpha.boot.trees.binned.morphospace, function(x) dispRity::dispRity(x, metric = c(median, centroids)))
Ornithuromorpha.boot.trees2.disparity.moc <- lapply(Ornithuromorpha.boot.trees2.binned.morphospace, function(x) dispRity::dispRity(x, metric = c(median, centroids)))
Ornithuromorpha.boot.sc.tree.disparity.moc <- dispRity::dispRity(Ornithuromorpha.boot.sc.tree.binned.morphospace, metric = c(median, centroids))
Ornithuromorpha.boot.sc.tree2.disparity.moc <- dispRity::dispRity(Ornithuromorpha.boot.sc.tree2.binned.morphospace, metric = c(median, centroids))
Ornithuromorpha.boot.tip.dated.tree.disparity.moc <- dispRity::dispRity(Ornithuromorpha.boot.tip.dated.tree.binned.morphospace, metric = c(median, centroids))

Plot sum of variance.

par(mfrow = c(2, 1))
plot(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = c(0, 0, 0), xlim = c(150, 66), ylim = c(0, 4.5), type = "n", xlab = "Time (Ma)", ylab = "Sum Of Variance", main = "Enantiornithes")
trees1 <- lapply(Enantiornithes.boot.trees.disparity.sov, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.1), lwd = 0.3)))
trees2 <- lapply(Enantiornithes.boot.trees2.disparity.sov, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.1), lwd = 0.3)))
trees3 <- apply(do.call(rbind, lapply(Enantiornithes.boot.sc.tree.disparity.sov$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.2), lwd = 0.3))
trees4 <- apply(do.call(rbind, lapply(Enantiornithes.boot.sc.tree2.disparity.sov$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.2), lwd = 0.3))
trees5 <- apply(do.call(rbind, lapply(Enantiornithes.boot.tip.dated.tree.disparity.sov$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Enantiornithes.boot.trees.disparity.sov, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Enantiornithes.boot.trees.disparity.sov, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "red", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Enantiornithes.boot.trees2.disparity.sov, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Enantiornithes.boot.trees2.disparity.sov, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "blue", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.sc.tree.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.sc.tree.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "red", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.sc.tree2.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.sc.tree2.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "blue", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.tip.dated.tree.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.tip.dated.tree.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)
plot(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = c(0, 0, 0), xlim = c(150, 66), ylim = c(0, 4.5), type = "n", xlab = "Time (Ma)", ylab = "Sum Of Variance", main = "Ornithuromorpha")
trees1 <- lapply(Ornithuromorpha.boot.trees.disparity.sov, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.1), lwd = 0.3)))
trees2 <- lapply(Ornithuromorpha.boot.trees2.disparity.sov, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.1), lwd = 0.3)))
trees3 <- apply(do.call(rbind, lapply(Ornithuromorpha.boot.sc.tree.disparity.sov$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.2), lwd = 0.3))
trees4 <- apply(do.call(rbind, lapply(Ornithuromorpha.boot.sc.tree2.disparity.sov$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.2), lwd = 0.3))
trees5 <- apply(do.call(rbind, lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.sov$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Ornithuromorpha.boot.trees.disparity.sov, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Ornithuromorpha.boot.trees.disparity.sov, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "red", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Ornithuromorpha.boot.trees2.disparity.sov, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Ornithuromorpha.boot.trees2.disparity.sov, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "blue", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.sc.tree.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.sc.tree.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "red", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.sc.tree2.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.sc.tree2.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "blue", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)

Plot sum of ranges.

par(mfrow = c(2, 1))
plot(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = c(0, 0, 0), xlim = c(150, 66), ylim = c(0, 35), type = "n", xlab = "Time (Ma)", ylab = "Sum Of Ranges", main = "Enantiornithes")
trees1 <- lapply(Enantiornithes.boot.trees.disparity.sor, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.1), lwd = 0.3)))
trees2 <- lapply(Enantiornithes.boot.trees2.disparity.sor, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.1), lwd = 0.3)))
trees3 <- apply(do.call(rbind, lapply(Enantiornithes.boot.sc.tree.disparity.sor$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.2), lwd = 0.3))
trees4 <- apply(do.call(rbind, lapply(Enantiornithes.boot.sc.tree2.disparity.sor$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.2), lwd = 0.3))
trees5 <- apply(do.call(rbind, lapply(Enantiornithes.boot.tip.dated.tree.disparity.sor$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Enantiornithes.boot.trees.disparity.sor, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Enantiornithes.boot.trees.disparity.sor, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "red", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Enantiornithes.boot.trees2.disparity.sor, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Enantiornithes.boot.trees2.disparity.sor, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "blue", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.sc.tree.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.sc.tree.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "red", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.sc.tree2.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.sc.tree2.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "blue", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.tip.dated.tree.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.tip.dated.tree.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)
plot(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = c(0, 0, 0), xlim = c(150, 66), ylim = c(0, 35), type = "n", xlab = "Time (Ma)", ylab = "Sum Of Ranges", main = "Ornithuromorpha")
trees1 <- lapply(Ornithuromorpha.boot.trees.disparity.sor, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.1), lwd = 0.3)))
trees2 <- lapply(Ornithuromorpha.boot.trees2.disparity.sor, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.1), lwd = 0.3)))
trees3 <- apply(do.call(rbind, lapply(Ornithuromorpha.boot.sc.tree.disparity.sor$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.2), lwd = 0.3))
trees4 <- apply(do.call(rbind, lapply(Ornithuromorpha.boot.sc.tree2.disparity.sor$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.2), lwd = 0.3))
trees5 <- apply(do.call(rbind, lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.sor$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Ornithuromorpha.boot.trees.disparity.sor, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Ornithuromorpha.boot.trees.disparity.sor, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "red", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Ornithuromorpha.boot.trees2.disparity.sor, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Ornithuromorpha.boot.trees2.disparity.sor, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "blue", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.sc.tree.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.sc.tree.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "red", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.sc.tree2.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.sc.tree2.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "blue", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)

Plot median centroid.

par(mfrow = c(2, 1))
plot(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = c(0, 0, 0), xlim = c(150, 66), ylim = c(0, 2), type = "n", xlab = "Time (Ma)", ylab = "Centroid Median", main = "Enantiornithes")
trees1 <- lapply(Enantiornithes.boot.trees.disparity.moc, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.1), lwd = 0.3)))
trees2 <- lapply(Enantiornithes.boot.trees2.disparity.moc, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.1), lwd = 0.3)))
trees3 <- apply(do.call(rbind, lapply(Enantiornithes.boot.sc.tree.disparity.moc$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.2), lwd = 0.3))
trees4 <- apply(do.call(rbind, lapply(Enantiornithes.boot.sc.tree2.disparity.moc$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.2), lwd = 0.3))
trees5 <- apply(do.call(rbind, lapply(Enantiornithes.boot.tip.dated.tree.disparity.moc$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Enantiornithes.boot.trees.disparity.moc, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Enantiornithes.boot.trees.disparity.moc, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "red", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Enantiornithes.boot.trees2.disparity.moc, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Enantiornithes.boot.trees2.disparity.moc, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "blue", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.sc.tree.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.sc.tree.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "red", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.sc.tree2.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.sc.tree2.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "blue", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.tip.dated.tree.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.tip.dated.tree.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)
plot(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = c(0, 0, 0), xlim = c(150, 66), ylim = c(0, 2), type = "n", xlab = "Time (Ma)", ylab = "Centroid Median", main = "Ornithuromorpha")
trees1 <- lapply(Ornithuromorpha.boot.trees.disparity.moc, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.1), lwd = 0.3)))
trees2 <- lapply(Ornithuromorpha.boot.trees2.disparity.moc, function(x) apply(do.call(rbind, lapply(x$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.1), lwd = 0.3)))
trees3 <- apply(do.call(rbind, lapply(Ornithuromorpha.boot.sc.tree.disparity.moc$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0, blue = 0, alpha = 0.2), lwd = 0.3))
trees4 <- apply(do.call(rbind, lapply(Ornithuromorpha.boot.sc.tree2.disparity.moc$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 0, green = 0, blue = 1, alpha = 0.2), lwd = 0.3))
trees5 <- apply(do.call(rbind, lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.moc$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Ornithuromorpha.boot.trees.disparity.moc, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Ornithuromorpha.boot.trees.disparity.moc, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "red", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Ornithuromorpha.boot.trees2.disparity.moc, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = apply(do.call(rbind, lapply(Ornithuromorpha.boot.trees2.disparity.moc, function(x) unlist(lapply(x$disparity, function(y) mean(unlist(y)))))), 2, mean), type = "l", col = "blue", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.sc.tree.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.sc.tree.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "red", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.sc.tree2.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.sc.tree2.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "blue", lwd = 8)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)

We can also do some PERMANOVA tests between our groups. First we will define a custom adonis function.

pairwise.adonis <- function(x, factors, sim.method, p.adjust.m) {
  co <- as.matrix(combn(unique(factors), 2))
  pairs <- c()
  F.Model <- c()
  R2 <- c()
  p.value <- c()
  for(elem in 1:ncol(co)) {
    ad <- vegan::adonis(x[factors %in% c(as.character(co[1, elem]), as.character(co[2, elem])), ] ~ factors[factors %in% c(as.character(co[1, elem]), as.character(co[2, elem]))], method = sim.method)
    pairs <- c(pairs,paste(co[1, elem], '&', co[2, elem]))
    F.Model <- c(F.Model, ad$aov.tab[1, 4])
    R2 <- c(R2, ad$aov.tab[1, 5])
    p.value <- c(p.value, ad$aov.tab[1, 6])
  }
  p.adjusted <- p.adjust(p.value, method = p.adjust.m)
  data.frame(pairs, F.Model, R2, p.value, p.adjusted)
}

We first need to create a factors vector for our taxon groups.

groups_as_factors <- as.factor(do.call(rbind, lapply(as.list(names(taxon_groups)), function(x) cbind(x, taxon_groups[[x]])))[, 1])
names(groups_as_factors) <- unname(unlist(taxon_groups))
groups_as_factors <- groups_as_factors[rownames(pcoa.data)]

And now we can do our pairwise adonis tests.

pairwise.adonis(pcoa.data, groups_as_factors, sim.method = "mahalanobis", p.adjust.m = "bonferroni")

Adding ancestral state estimates to the data.

We can repeat many of these tests but including ancestral state estimates in the sample. We can start with Claddis.

First we must remake our trees by only excluding crown taxa.

stem.time.trees <- lapply(time.trees, function(x) Claddis::drop_time_tip(x, tip_names = crown.taxa.to.remove))
stem.time.trees2 <- lapply(time.trees, function(x) Claddis::drop_time_tip(x, tip_names = crown.taxa.to.remove))
stem.tip.dated.tree <- Claddis::drop_time_tip(tip.dated.tree, tip_names = crown.taxa.to.remove)

Now we can do a phylogenetic ordination where discrete ancestral character estimations are made for internal nodes using our tip-dated tree.

stem.tip.dated.tree.pcoa.data <- Claddis::ordinate_cladistic_matrix(stem.nexus.data, time_tree = stem.tip.dated.tree)

We can make simple plots of the resulting space.

Claddis::plot_morphospace(stem.tip.dated.tree.pcoa.data, x_axis = 1, y_axis = 2, plot_taxon_names = TRUE, plot_convex_hulls = TRUE, taxon_groups = taxon_groups)

## [1] "Warning: The following taxa were removed from taxon_groups as they do not appear in pcoa_input: Shangyang. You may wish to double check this makes sense (e.g., because of incomplete taxa being removed by trim_matrix) and is not due to a typographical or other error which means names are not an exact match."

Claddis::plot_morphospace(stem.tip.dated.tree.pcoa.data, x_axis = 1, y_axis = 3, plot_taxon_names = TRUE, plot_convex_hulls = TRUE, taxon_groups = taxon_groups)

## [1] "Warning: The following taxa were removed from taxon_groups as they do not appear in pcoa_input: Shangyang. You may wish to double check this makes sense (e.g., because of incomplete taxa being removed by trim_matrix) and is not due to a typographical or other error which means names are not an exact match."

Claddis::plot_morphospace(stem.tip.dated.tree.pcoa.data, x_axis = 2, y_axis = 3, plot_taxon_names = TRUE, plot_convex_hulls = TRUE, taxon_groups = taxon_groups)

## [1] "Warning: The following taxa were removed from taxon_groups as they do not appear in pcoa_input: Shangyang. You may wish to double check this makes sense (e.g., because of incomplete taxa being removed by trim_matrix) and is not due to a typographical or other error which means names are not an exact match."

Or multiple plots.

Claddis::plot_multi_morphospace(stem.tip.dated.tree.pcoa.data, n_axes = 8, plot_convex_hulls = TRUE, taxon_groups = taxon_groups)

And calculate sum of ranges…

c(Enantiornithes = sum(abs(apply(apply(stem.tip.dated.tree.pcoa.data$vectors[intersect(rownames(stem.tip.dated.tree.pcoa.data$vectors), Enantiornithes), ], 2, range), 2, diff))), Ornithuromorpha = sum(abs(apply(apply(stem.tip.dated.tree.pcoa.data$vectors[intersect(rownames(stem.tip.dated.tree.pcoa.data$vectors), Ornithuromorpha), ], 2, range), 2, diff))), NonOrnithothoraces = sum(abs(apply(apply(stem.tip.dated.tree.pcoa.data$vectors[intersect(rownames(stem.tip.dated.tree.pcoa.data$vectors), NonOrnithothoraces), ], 2, range), 2, diff))))

…sum of variances…

c(Enantiornithes = sum(apply(stem.pcoa.data$vectors[intersect(rownames(stem.tip.dated.tree.pcoa.data$vectors), Enantiornithes), ], 2, var)), Ornithuromorpha = sum(apply(stem.pcoa.data$vectors[intersect(rownames(stem.tip.dated.tree.pcoa.data$vectors), Ornithuromorpha), ], 2, var)), NonOrnithothoraces = sum(apply(stem.pcoa.data$vectors[intersect(rownames(stem.tip.dated.tree.pcoa.data$vectors), NonOrnithothoraces), ], 2, var)))

…product of ranges….

c(Enantiornithes = prod(abs(apply(apply(stem.pcoa.data$vectors[intersect(rownames(stem.tip.dated.tree.pcoa.data$vectors), Enantiornithes), ], 2, range), 2, diff))), Ornithuromorpha = prod(abs(apply(apply(stem.pcoa.data$vectors[intersect(rownames(stem.tip.dated.tree.pcoa.data$vectors), Ornithuromorpha), ], 2, range), 2, diff))), NonOrnithothoraces = prod(abs(apply(apply(stem.pcoa.data$vectors[intersect(rownames(stem.tip.dated.tree.pcoa.data$vectors), NonOrnithothoraces), ], 2, range), 2, diff))))

…or product of variances.

c(Enantiornithes = prod(apply(stem.pcoa.data$vectors[intersect(rownames(stem.tip.dated.tree.pcoa.data$vectors), Enantiornithes), ], 2, var)), Ornithuromorpha = prod(apply(stem.pcoa.data$vectors[intersect(rownames(stem.tip.dated.tree.pcoa.data$vectors), Ornithuromorpha), ], 2, var)), NonOrnithothoraces = prod(apply(stem.pcoa.data$vectors[intersect(rownames(stem.tip.dated.tree.pcoa.data$vectors), NonOrnithothoraces), ], 2, var)))

Next we can repeat our dispRity analyses. First we need to trim down the tree again.

stem.tip.dated.tree <- stem.tip.dated.tree.pcoa.data$time_tree

We also need to add node labels to the tree.

stem.tip.dated.tree$node.label <- paste0("n", 1:stem.tip.dated.tree$Nnode)

And we can create a new phylogenetic taxon groups object too.

phylo_taxon_groups <- lapply(taxon_groups, function(x) c(stem.tip.dated.tree$tip.label[intersect(1:ape::Ntip(stem.tip.dated.tree), stem.tip.dated.tree$edge[Claddis::find_descendant_edges(n = Claddis::find_mrca(intersect(x, stem.tip.dated.tree$tip.label), stem.tip.dated.tree), tree = stem.tip.dated.tree), 2])], stem.tip.dated.tree$node.label[unique(stem.tip.dated.tree$edge[Claddis::find_descendant_edges(n = Claddis::find_mrca(intersect(x, stem.tip.dated.tree$tip.label), stem.tip.dated.tree), tree = stem.tip.dated.tree), 1] - ape::Ntip(stem.tip.dated.tree))]))
phylo_taxon_groups$NonOrnithothoraces <- setdiff(c(stem.tip.dated.tree$tip.label, stem.tip.dated.tree$node.label), unlist(phylo_taxon_groups[c("Enantiornithes", "Ornithuromorpha")]))

And we will only use the first 40 axes.

stem.tip.dated.tree.pcoa.data <- as.data.frame(stem.tip.dated.tree.pcoa.data$vectors[, 1:40])
stem.tip.dated.tree.pcoa.data <- data.matrix(stem.tip.dated.tree.pcoa.data)
rownames(stem.tip.dated.tree.pcoa.data)[match(as.character((ape::Ntip(stem.tip.dated.tree) + 1):(ape::Ntip(stem.tip.dated.tree) + stem.tip.dated.tree$Nnode)), rownames(stem.tip.dated.tree.pcoa.data))] <- paste0("n", 1:stem.tip.dated.tree$Nnode)

And as we are going to use the dispRity package convert our ages to the preferred format.

disprity.ages.data <- tip.ages[stem.tip.dated.tree$tip.label, ]
colnames(disprity.ages.data) <- c("FAD", "LAD")

And set some new time bins for the time series analyses.

time.bins <- c(170, 145, 129.4, 100.5, 86.3, 66)

We can now use our various phylogenetc hypotheses and time bin scheme to create dispRity objects.

stem.tip.dated.tree.binned.morphospace <- dispRity::chrono.subsets(data = stem.tip.dated.tree.pcoa.data, tree = stem.tip.dated.tree, method = "discrete", time = time.bins, inc.nodes = TRUE, FADLAD = disprity.ages.data)

We can also use dispRity to bootstrap our data and give a better sense of our uncertianty.

stem.tip.dated.tree.boot.bin.morphospace <- dispRity::boot.matrix(stem.tip.dated.tree.binned.morphospace, bootstrap = 1000)

Now we can calculate some disparity metrics again, such as sum of variances.

stem.tip.dated.tree.boot.bin.morphospace.sov <- dispRity::dispRity(stem.tip.dated.tree.boot.bin.morphospace, metric = c(sum, variances))

And visualise them.

plot(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = c(0, 0, 0, 0, 0), xlim = c(170, 66), ylim = c(0, 12), type = "n", xlab = "Time (Ma)", ylab = "Sum Of Variance")
trees5 <- apply(do.call(rbind, lapply(stem.tip.dated.tree.boot.bin.morphospace.sov$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.tip.dated.tree.boot.bin.morphospace.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.tip.dated.tree.boot.bin.morphospace.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)

Or sum of ranges.

stem.tip.dated.tree.boot.bin.morphospace.sor <- dispRity::dispRity(stem.tip.dated.tree.boot.bin.morphospace, metric = c(sum, ranges))

And visualise them.

plot(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = c(0, 0, 0, 0, 0), xlim = c(170, 66), ylim = c(0, 100), type = "n", xlab = "Time (Ma)", ylab = "Sum Of Ranges")
trees5 <- apply(do.call(rbind, lapply(stem.tip.dated.tree.boot.bin.morphospace.sor$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.tip.dated.tree.boot.bin.morphospace.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.tip.dated.tree.boot.bin.morphospace.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)

Or centroid medians.

stem.tip.dated.tree.boot.bin.morphospace.moc <- dispRity::dispRity(stem.tip.dated.tree.boot.bin.morphospace, metric = c(median, centroids))

And visualise them.

plot(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = c(0, 0, 0, 0, 0), xlim = c(170, 66), ylim = c(0, 3.5), type = "n", xlab = "Time (Ma)", ylab = "Centroid Median")
trees5 <- apply(do.call(rbind, lapply(stem.tip.dated.tree.boot.bin.morphospace.moc$disparity, function(y) unlist(y))), 2, function(z) points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.tip.dated.tree.boot.bin.morphospace.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = time.bins[1:(length(time.bins) - 1)] + diff(time.bins) / 2, y = unlist(lapply(stem.tip.dated.tree.boot.bin.morphospace.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)

We can also subset by group.

stem.pcoa.data.groups <- dispRity::custom.subsets(stem.tip.dated.tree.pcoa.data, group = lapply(phylo_taxon_groups, function(x) intersect(c(stem.tip.dated.tree$tip.label, stem.tip.dated.tree$node.label), x)))

Bootstrap.

stem.pcoa.data.groups.bootstrap <- dispRity::boot.matrix(stem.pcoa.data.groups, bootstrap = 1000)

And calculate our disparity metrics.

stem.pcoa.data.groups.bootstrap.sov <- dispRity::dispRity(stem.pcoa.data.groups.bootstrap, metric = c(sum, variances))
stem.pcoa.data.groups.bootstrap.sor <- dispRity::dispRity(stem.pcoa.data.groups.bootstrap, metric = c(sum, ranges))
stem.pcoa.data.groups.bootstrap.moc <- dispRity::dispRity(stem.pcoa.data.groups.bootstrap, metric = c(median, centroids))

And plot the results.

par(mfrow=c(1,3))
boxplot(lapply(stem.pcoa.data.groups.bootstrap.sov$disparity, unlist), ylab = "Sum of Variance", main = "Sum of Variance")
boxplot(lapply(stem.pcoa.data.groups.bootstrap.sor$disparity, unlist), ylab = "Sum of Ranges", main = "Sum of Ranges")
boxplot(lapply(stem.pcoa.data.groups.bootstrap.moc$disparity, unlist), ylab = "Centroid Median", main = "Centroid Median")

We can also repeat the Enantiornithes and Ornithuromorpha time series dispRity analyses from above, but this time including ancestral state estimates in the sample.

First we need to remake the dispRity objects.

Enantiornithes.tip.dated.tree.binned.morphospace <- dispRity::chrono.subsets(data = stem.tip.dated.tree.pcoa.data[c(Claddis::drop_time_tip(stem.tip.dated.tree, setdiff(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Enantiornithes))$tip.label, Claddis::drop_time_tip(stem.tip.dated.tree, setdiff(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Enantiornithes))$node.label), ], tree = Claddis::drop_time_tip(stem.tip.dated.tree, setdiff(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Enantiornithes)), method = "discrete", time = comparison_time_bins, inc.nodes = TRUE, FADLAD = disprity.ages.data[intersect(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Enantiornithes), ])
Ornithuromorpha.tip.dated.tree.binned.morphospace <- dispRity::chrono.subsets(data = stem.tip.dated.tree.pcoa.data[c(Claddis::drop_time_tip(stem.tip.dated.tree, setdiff(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Ornithuromorpha))$tip.label, Claddis::drop_time_tip(stem.tip.dated.tree, setdiff(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Ornithuromorpha))$node.label), ], tree = Claddis::drop_time_tip(stem.tip.dated.tree, setdiff(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Ornithuromorpha)), method = "discrete", time = comparison_time_bins, inc.nodes = TRUE, FADLAD = disprity.ages.data[intersect(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Ornithuromorpha), ])

And bootstrap the data.

Enantiornithes.boot.tip.dated.tree.binned.morphospace <- dispRity::boot.matrix(Enantiornithes.tip.dated.tree.binned.morphospace, bootstrap = 1000)
Ornithuromorpha.boot.tip.dated.tree.binned.morphospace <- dispRity::boot.matrix(Ornithuromorpha.tip.dated.tree.binned.morphospace, bootstrap = 1000)

Calculate sum of variances.

Enantiornithes.boot.tip.dated.tree.disparity.sov <- dispRity::dispRity(Enantiornithes.boot.tip.dated.tree.binned.morphospace, metric = c(sum, variances))
Ornithuromorpha.boot.tip.dated.tree.disparity.sov <- dispRity::dispRity(Ornithuromorpha.boot.tip.dated.tree.binned.morphospace, metric = c(sum, variances))

Calculate sum of ranges.

Enantiornithes.boot.tip.dated.tree.disparity.sor <- dispRity::dispRity(Enantiornithes.boot.tip.dated.tree.binned.morphospace, metric = c(sum, ranges))
Ornithuromorpha.boot.tip.dated.tree.disparity.sor <- dispRity::dispRity(Ornithuromorpha.boot.tip.dated.tree.binned.morphospace, metric = c(sum, ranges))

Calculate median centroids.

Enantiornithes.boot.tip.dated.tree.disparity.moc <- dispRity::dispRity(Enantiornithes.boot.tip.dated.tree.binned.morphospace, metric = c(median, centroids))
Ornithuromorpha.boot.tip.dated.tree.disparity.moc <- dispRity::dispRity(Ornithuromorpha.boot.tip.dated.tree.binned.morphospace, metric = c(median, centroids))

Plot sum of variance.

par(mfrow = c(2, 1))
plot(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = c(0, 0, 0), xlim = c(150, 66), ylim = c(0, 10), type = "n", xlab = "Time (Ma)", ylab = "Sum Of Variance", main = "Enantiornithes")
trees5 <- apply(do.call(rbind, lapply(Enantiornithes.boot.tip.dated.tree.disparity.sov$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.tip.dated.tree.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.tip.dated.tree.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)
plot(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = c(0, 0, 0), xlim = c(150, 66), ylim = c(0, 10), type = "n", xlab = "Time (Ma)", ylab = "Sum Of Variance", main = "Ornithuromorpha")
trees5 <- apply(do.call(rbind, lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.sov$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.sov$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)

Plot sum of ranges.

par(mfrow = c(2, 1))
plot(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = c(0, 0, 0), xlim = c(150, 66), ylim = c(0, 75), type = "n", xlab = "Time (Ma)", ylab = "Sum Of Ranges", main = "Enantiornithes")
trees5 <- apply(do.call(rbind, lapply(Enantiornithes.boot.tip.dated.tree.disparity.sor$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.tip.dated.tree.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.tip.dated.tree.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)
plot(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = c(0, 0, 0), xlim = c(150, 66), ylim = c(0, 75), type = "n", xlab = "Time (Ma)", ylab = "Sum Of Ranges", main = "Ornithuromorpha")
trees5 <- apply(do.call(rbind, lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.sor$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.sor$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)

Plot median centroid.

par(mfrow = c(2, 1))
plot(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = c(0, 0, 0), xlim = c(150, 66), ylim = c(0, 3), type = "n", xlab = "Time (Ma)", ylab = "Centroid Median", main = "Enantiornithes")
trees5 <- apply(do.call(rbind, lapply(Enantiornithes.boot.tip.dated.tree.disparity.moc$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.tip.dated.tree.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Enantiornithes.boot.tip.dated.tree.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)
plot(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = c(0, 0, 0), xlim = c(150, 66), ylim = c(0, 3), type = "n", xlab = "Time (Ma)", ylab = "Centroid Median", main = "Ornithuromorpha")
trees5 <- apply(do.call(rbind, lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.moc$disparity, function(y) unlist(y))), 2, function(z) points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = z, type = "l", col = rgb(red = 1, green = 0.25, blue = 0, alpha = 0.2), lwd = 0.3))
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "black", lwd = 10)
points(x = comparison_time_bins[1:(length(comparison_time_bins) - 1)] + diff(comparison_time_bins) / 2, y = unlist(lapply(Ornithuromorpha.boot.tip.dated.tree.disparity.moc$disparity, function(x) mean(unlist(x)))), type = "l", col = "orange", lwd = 8)

We can also repeat the PERMANOVA test before. First we can create new taxon groups.

taxon_groups_inc_ancestors <- list(NonOrnithothoracesIncAncestors = setdiff(rownames(stem.tip.dated.tree.pcoa.data), c(Claddis::drop_time_tip(stem.tip.dated.tree, setdiff(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Enantiornithes))$tip.label, Claddis::drop_time_tip(stem.tip.dated.tree, setdiff(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Enantiornithes))$node.label, Claddis::drop_time_tip(stem.tip.dated.tree, setdiff(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Ornithuromorpha))$tip.label, Claddis::drop_time_tip(stem.tip.dated.tree, setdiff(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Ornithuromorpha))$node.label)), EnantiornithesIncAncestors = c(Claddis::drop_time_tip(stem.tip.dated.tree, setdiff(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Enantiornithes))$tip.label, Claddis::drop_time_tip(stem.tip.dated.tree, setdiff(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Enantiornithes))$node.label), OrnithuromorphaIncAncestors = c(Claddis::drop_time_tip(stem.tip.dated.tree, setdiff(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Ornithuromorpha))$tip.label, Claddis::drop_time_tip(stem.tip.dated.tree, setdiff(stem.tip.dated.tree$tip.label, phylo_taxon_groups$Ornithuromorpha))$node.label))
groups_inc_ancestors_as_factors <- as.factor(do.call(rbind, lapply(as.list(names(taxon_groups_inc_ancestors)), function(x) cbind(x, taxon_groups_inc_ancestors[[x]])))[, 1])
names(groups_inc_ancestors_as_factors) <- unname(unlist(taxon_groups_inc_ancestors))
groups_inc_ancestors_as_factors <- groups_inc_ancestors_as_factors[rownames(stem.tip.dated.tree.pcoa.data)]

Perform test.

pairwise.adonis(stem.tip.dated.tree.pcoa.data, groups_inc_ancestors_as_factors, sim.method = "mahalanobis", p.adjust.m = "bonferroni")

Reanalysing by skeletal subregion

We can also examine individual subregions of the skeleton.

First we can define these as the character sets they contain.

skull.characters <- c(1:48, 262, 263, 267, 268, 275)
axial.characters <- c(49:81, 269, 270)
pectoral.characters <- c(82:120, 246:249, 251, 252, 271, 274, 276)
forelimb.characters <- c(121:177, 250, 253, 265, 277:280)
pelvic.characters <- c(178:198, 272, 273)
hindlimb.characters <- c(199:243, 254:261, 264, 266)

And use these to make pruned versions of our cladistic matrix.

skull.stem.nexus.data <- Claddis::prune_cladistic_matrix(stem.nexus.data, characters2prune = setdiff(1:ncol(stem.nexus.data$matrix_1$matrix), skull.characters))
axial.stem.nexus.data <- Claddis::prune_cladistic_matrix(stem.nexus.data, characters2prune = setdiff(1:ncol(stem.nexus.data$matrix_1$matrix), axial.characters))
pectoral.stem.nexus.data <- Claddis::prune_cladistic_matrix(stem.nexus.data, characters2prune = setdiff(1:ncol(stem.nexus.data$matrix_1$matrix), pectoral.characters))
forelimb.stem.nexus.data <- Claddis::prune_cladistic_matrix(stem.nexus.data, characters2prune = setdiff(1:ncol(stem.nexus.data$matrix_1$matrix), forelimb.characters))
pelvic.stem.nexus.data <- Claddis::prune_cladistic_matrix(stem.nexus.data, characters2prune = setdiff(1:ncol(stem.nexus.data$matrix_1$matrix), pelvic.characters))
hindlimb.stem.nexus.data <- Claddis::prune_cladistic_matrix(stem.nexus.data, characters2prune = setdiff(1:ncol(stem.nexus.data$matrix_1$matrix), hindlimb.characters))

And we can now ordinate the data.

skull.pcoa.data <- Claddis::ordinate_cladistic_matrix(skull.stem.nexus.data, distance_metric = "mord")
axial.pcoa.data <- Claddis::ordinate_cladistic_matrix(axial.stem.nexus.data, distance_metric = "mord")
pectoral.pcoa.data <- Claddis::ordinate_cladistic_matrix(pectoral.stem.nexus.data, distance_metric = "mord")
forelimb.pcoa.data <- Claddis::ordinate_cladistic_matrix(forelimb.stem.nexus.data, distance_metric = "mord")
pelvic.pcoa.data <- Claddis::ordinate_cladistic_matrix(pelvic.stem.nexus.data, distance_metric = "mord")
hindlimb.pcoa.data <- Claddis::ordinate_cladistic_matrix(hindlimb.stem.nexus.data, distance_metric = "mord")

(Note that many more taxa have to be removed for each subregion than for the total data set.)

We can now convert thesd into dispRity objects, subsettng by time.

skull.chrono.dispRity <- dispRity::chrono.subsets(data = skull.pcoa.data$vectors, tree = Claddis::drop_time_tip(tip.dated.tree, setdiff(tip.dated.tree$tip.label, rownames(skull.pcoa.data$vectors))), method = "discrete", time = time.bins, inc.nodes = FALSE, FADLAD = tip.ages[rownames(skull.pcoa.data$vectors), ])
axial.chrono.dispRity <- dispRity::chrono.subsets(data = axial.pcoa.data$vectors, tree = Claddis::drop_time_tip(tip.dated.tree, setdiff(tip.dated.tree$tip.label, rownames(axial.pcoa.data$vectors))), method = "discrete", time = time.bins, inc.nodes = FALSE, FADLAD = tip.ages[rownames(axial.pcoa.data$vectors), ])
pectoral.chrono.dispRity <- dispRity::chrono.subsets(data = pectoral.pcoa.data$vectors, tree = Claddis::drop_time_tip(tip.dated.tree, setdiff(tip.dated.tree$tip.label, rownames(pectoral.pcoa.data$vectors))), method = "discrete", time = time.bins, inc.nodes = FALSE, FADLAD = tip.ages[rownames(pectoral.pcoa.data$vectors), ])
forelimb.chrono.dispRity <- dispRity::chrono.subsets(data = forelimb.pcoa.data$vectors, tree = Claddis::drop_time_tip(tip.dated.tree, setdiff(tip.dated.tree$tip.label, rownames(forelimb.pcoa.data$vectors))), method = "discrete", time = time.bins, inc.nodes = FALSE, FADLAD = tip.ages[rownames(forelimb.pcoa.data$vectors), ])
pelvic.chrono.dispRity <- dispRity::chrono.subsets(data = pelvic.pcoa.data$vectors, tree = Claddis::drop_time_tip(tip.dated.tree, setdiff(tip.dated.tree$tip.label, rownames(pelvic.pcoa.data$vectors))), method = "discrete", time = time.bins, inc.nodes = FALSE, FADLAD = tip.ages[rownames(pelvic.pcoa.data$vectors), ])
hindlimb.chrono.dispRity <- dispRity::chrono.subsets(data = hindlimb.pcoa.data$vectors, tree = Claddis::drop_time_tip(tip.dated.tree, setdiff(tip.dated.tree$tip.label, rownames(hindlimb.pcoa.data$vectors))), method = "discrete", time = time.bins, inc.nodes = FALSE, FADLAD = tip.ages[rownames(hindlimb.pcoa.data$vectors), ])

Or by group.

skull.group.dispRity <- dispRity::custom.subsets(skull.pcoa.data$vectors, group = list(NonOrnithothoraces = match(intersect(NonOrnithothoraces, rownames(skull.pcoa.data$vectors)), rownames(skull.pcoa.data$vectors)), Enantiornithes = match(intersect(Enantiornithes, rownames(skull.pcoa.data$vectors)), rownames(skull.pcoa.data$vectors)), Ornithuromorpha = match(intersect(Ornithuromorpha, rownames(skull.pcoa.data$vectors)), rownames(skull.pcoa.data$vectors))))
axial.group.dispRity <- dispRity::custom.subsets(axial.pcoa.data$vectors, group = list(NonOrnithothoraces = match(intersect(NonOrnithothoraces, rownames(axial.pcoa.data$vectors)), rownames(axial.pcoa.data$vectors)), Enantiornithes = match(intersect(Enantiornithes, rownames(axial.pcoa.data$vectors)), rownames(axial.pcoa.data$vectors)), Ornithuromorpha = match(intersect(Ornithuromorpha, rownames(axial.pcoa.data$vectors)), rownames(axial.pcoa.data$vectors))))
pectoral.group.dispRity <- dispRity::custom.subsets(pectoral.pcoa.data$vectors, group = list(NonOrnithothoraces = match(intersect(NonOrnithothoraces, rownames(pectoral.pcoa.data$vectors)), rownames(pectoral.pcoa.data$vectors)), Enantiornithes = match(intersect(Enantiornithes, rownames(pectoral.pcoa.data$vectors)), rownames(pectoral.pcoa.data$vectors)), Ornithuromorpha = match(intersect(Ornithuromorpha, rownames(pectoral.pcoa.data$vectors)), rownames(pectoral.pcoa.data$vectors))))
forelimb.group.dispRity <- dispRity::custom.subsets(forelimb.pcoa.data$vectors, group = list(NonOrnithothoraces = match(intersect(NonOrnithothoraces, rownames(forelimb.pcoa.data$vectors)), rownames(forelimb.pcoa.data$vectors)), Enantiornithes = match(intersect(Enantiornithes, rownames(forelimb.pcoa.data$vectors)), rownames(forelimb.pcoa.data$vectors)), Ornithuromorpha = match(intersect(Ornithuromorpha, rownames(forelimb.pcoa.data$vectors)), rownames(forelimb.pcoa.data$vectors))))
pelvic.group.dispRity <- dispRity::custom.subsets(pelvic.pcoa.data$vectors, group = list(NonOrnithothoraces = match(intersect(NonOrnithothoraces, rownames(pelvic.pcoa.data$vectors)), rownames(pelvic.pcoa.data$vectors)), Enantiornithes = match(intersect(Enantiornithes, rownames(pelvic.pcoa.data$vectors)), rownames(pelvic.pcoa.data$vectors)), Ornithuromorpha = match(intersect(Ornithuromorpha, rownames(pelvic.pcoa.data$vectors)), rownames(pelvic.pcoa.data$vectors))))
hindlimb.group.dispRity <- dispRity::custom.subsets(hindlimb.pcoa.data$vectors, group = list(NonOrnithothoraces = match(intersect(NonOrnithothoraces, rownames(hindlimb.pcoa.data$vectors)), rownames(hindlimb.pcoa.data$vectors)), Enantiornithes = match(intersect(Enantiornithes, rownames(hindlimb.pcoa.data$vectors)), rownames(hindlimb.pcoa.data$vectors)), Ornithuromorpha = match(intersect(Ornithuromorpha, rownames(hindlimb.pcoa.data$vectors)), rownames(hindlimb.pcoa.data$vectors))))

And bootstrap the data.

skull.chrono.dispRity.booted <- dispRity::boot.matrix(skull.chrono.dispRity, bootstrap = 1000)
axial.chrono.dispRity.booted <- dispRity::boot.matrix(axial.chrono.dispRity, bootstrap = 1000)
pectoral.chrono.dispRity.booted <- dispRity::boot.matrix(pectoral.chrono.dispRity, bootstrap = 1000)
forelimb.chrono.dispRity.booted <- dispRity::boot.matrix(forelimb.chrono.dispRity, bootstrap = 1000)
pelvic.chrono.dispRity.booted <- dispRity::boot.matrix(pelvic.chrono.dispRity, bootstrap = 1000)
hindlimb.chrono.dispRity.booted <- dispRity::boot.matrix(hindlimb.chrono.dispRity, bootstrap = 1000)
skull.group.dispRity.booted <- dispRity::boot.matrix(skull.group.dispRity, bootstrap = 1000)
axial.group.dispRity.booted <- dispRity::boot.matrix(axial.group.dispRity, bootstrap = 1000)
pectoral.group.dispRity.booted <- dispRity::boot.matrix(pectoral.group.dispRity, bootstrap = 1000)
forelimb.group.dispRity.booted <- dispRity::boot.matrix(forelimb.group.dispRity, bootstrap = 1000)
pelvic.group.dispRity.booted <- dispRity::boot.matrix(pelvic.group.dispRity, bootstrap = 1000)
hindlimb.group.dispRity.booted <- dispRity::boot.matrix(hindlimb.group.dispRity, bootstrap = 1000)

We can now plot time series, of sum of variances.

plot(x = 0, y = 0, type = "n", xlim = c(170, 66), ylim = c(0, 25), xlab = "Time (Ma)", ylab = "Sum of Variances")
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(skull.chrono.dispRity.booted, metric = c(sum, variances))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[1])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(axial.chrono.dispRity.booted, metric = c(sum, variances))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[2])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(pectoral.chrono.dispRity.booted, metric = c(sum, variances))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[3])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(forelimb.chrono.dispRity.booted, metric = c(sum, variances))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[4])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(pelvic.chrono.dispRity.booted, metric = c(sum, variances))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[5])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(hindlimb.chrono.dispRity.booted, metric = c(sum, variances))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[6])
legend(x = 170, y = 25, legend = c("Skull", "Axial", "Pectoral", "Forelimb", "Pelvic", "Hindlimb"), fill = NULL, col = topo.colors(n = 6), border = "white", lty = 1, lwd = 4, pch = NULL, bg = "white")

Sum of ranges.

plot(x = 0, y = 0, type = "n", xlim = c(170, 66), ylim = c(0, 200), xlab = "Time (Ma)", ylab = "Sum of Ranges")
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(skull.chrono.dispRity.booted, metric = c(sum, ranges))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[1])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(axial.chrono.dispRity.booted, metric = c(sum, ranges))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[2])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(pectoral.chrono.dispRity.booted, metric = c(sum, ranges))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[3])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(forelimb.chrono.dispRity.booted, metric = c(sum, ranges))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[4])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(pelvic.chrono.dispRity.booted, metric = c(sum, ranges))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[5])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(hindlimb.chrono.dispRity.booted, metric = c(sum, ranges))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[6])
legend(x = 170, y = 200, legend = c("Skull", "Axial", "Pectoral", "Forelimb", "Pelvic", "Hindlimb"), fill = NULL, col = topo.colors(n = 6), border = "white", lty = 1, lwd = 4, pch = NULL, bg = "white")

Or centroid medians.

plot(x = 0, y = 0, type = "n", xlim = c(170, 66), ylim = c(0, 6), xlab = "Time (Ma)", ylab = "Median Centroid")
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(skull.chrono.dispRity.booted, metric = c(median, centroids))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[1])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(axial.chrono.dispRity.booted, metric = c(median, centroids))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[2])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(pectoral.chrono.dispRity.booted, metric = c(median, centroids))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[3])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(forelimb.chrono.dispRity.booted, metric = c(median, centroids))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[4])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(pelvic.chrono.dispRity.booted, metric = c(median, centroids))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[5])
points(x = time.bins[1:(length(time.bins) - 1)] + (diff(time.bins) / 2), y = unlist(lapply(dispRity::dispRity(hindlimb.chrono.dispRity.booted, metric = c(median, centroids))$disparity, function(x) mean(unlist(x)))), type = "l", lwd = 4, col = topo.colors(n = 6)[6])
legend(x = 170, y = 6, legend = c("Skull", "Axial", "Pectoral", "Forelimb", "Pelvic", "Hindlimb"), fill = NULL, col = topo.colors(n = 6), border = "white", lty = 1, lwd = 4, pch = NULL, bg = "white")

And now the same by groups, first sum of variances.

par(mfrow = c(2, 3))
boxplot(lapply(dispRity::dispRity(skull.group.dispRity.booted, metric = c(sum, variances))$disparity, function(x) unlist(x)), main = "Skull", ylab = "Sum of Variances")
boxplot(lapply(dispRity::dispRity(axial.group.dispRity.booted, metric = c(sum, variances))$disparity, function(x) unlist(x)), main = "Axial", ylab = "Sum of Variances")
boxplot(lapply(dispRity::dispRity(pectoral.group.dispRity.booted, metric = c(sum, variances))$disparity, function(x) unlist(x)), main = "Pectoral", ylab = "Sum of Variances")
boxplot(lapply(dispRity::dispRity(forelimb.group.dispRity.booted, metric = c(sum, variances))$disparity, function(x) unlist(x)), main = "Forelimb", ylab = "Sum of Variances")
boxplot(lapply(dispRity::dispRity(pelvic.group.dispRity.booted, metric = c(sum, variances))$disparity, function(x) unlist(x)), main = "Pelvic", ylab = "Sum of Variances")
boxplot(lapply(dispRity::dispRity(hindlimb.group.dispRity.booted, metric = c(sum, variances))$disparity, function(x) unlist(x)), main = "Hindlimb", ylab = "Sum of Variances")

Sum of ranges.

par(mfrow = c(2, 3))
boxplot(lapply(dispRity::dispRity(skull.group.dispRity.booted, metric = c(sum, ranges))$disparity, function(x) unlist(x)), main = "Skull", ylab = "Sum of Ranges")
boxplot(lapply(dispRity::dispRity(axial.group.dispRity.booted, metric = c(sum, ranges))$disparity, function(x) unlist(x)), main = "Axial", ylab = "Sum of Ranges")
boxplot(lapply(dispRity::dispRity(pectoral.group.dispRity.booted, metric = c(sum, ranges))$disparity, function(x) unlist(x)), main = "Pectoral", ylab = "Sum of Ranges")
boxplot(lapply(dispRity::dispRity(forelimb.group.dispRity.booted, metric = c(sum, ranges))$disparity, function(x) unlist(x)), main = "Forelimb", ylab = "Sum of Ranges")
boxplot(lapply(dispRity::dispRity(pelvic.group.dispRity.booted, metric = c(sum, ranges))$disparity, function(x) unlist(x)), main = "Pelvic", ylab = "Sum of Ranges")
boxplot(lapply(dispRity::dispRity(hindlimb.group.dispRity.booted, metric = c(sum, ranges))$disparity, function(x) unlist(x)), main = "Hindlimb", ylab = "Sum of Ranges")

Centroid medians.

par(mfrow = c(2, 3))
boxplot(lapply(dispRity::dispRity(skull.group.dispRity.booted, metric = c(median, centroids))$disparity, function(x) unlist(x)), main = "Skull", ylab = "Centroid Medians")
boxplot(lapply(dispRity::dispRity(axial.group.dispRity.booted, metric = c(median, centroids))$disparity, function(x) unlist(x)), main = "Axial", ylab = "Centroid Medians")
boxplot(lapply(dispRity::dispRity(pectoral.group.dispRity.booted, metric = c(median, centroids))$disparity, function(x) unlist(x)), main = "Pectoral", ylab = "Centroid Medians")
boxplot(lapply(dispRity::dispRity(forelimb.group.dispRity.booted, metric = c(median, centroids))$disparity, function(x) unlist(x)), main = "Forelimb", ylab = "Centroid Medians")
boxplot(lapply(dispRity::dispRity(pelvic.group.dispRity.booted, metric = c(median, centroids))$disparity, function(x) unlist(x)), main = "Pelvic", ylab = "Centroid Medians")
boxplot(lapply(dispRity::dispRity(hindlimb.group.dispRity.booted, metric = c(median, centroids))$disparity, function(x) unlist(x)), main = "Hindlimb", ylab = "Centroid Medians")

Effect of body mass on morphospace

Brougham and Campione Brougham and Campione (2020) showed that body size can strongly influence discrete character morphospaces in dinosaurs and we can use the same technique to explore whether a similar effect is seen in birds.

First we need to import body mass estimates for our taxa (NB: these are not available for all species).

mass <- matrix(data = c(143.47, 1084.33, 507.74, 825.25, 1982.61, 358.95, 138.4, 237.09, 173.41, 420.61, 112.49, 56.29, 81.42, 57.01, 88.85, 116.95, 43.01, 102.83, 524.68, 347.78, 111.39, 488.41, 188.18, 63.78, 22.69, 40.02, 38.3, 238.84, 281.34, 221.68, 220.01, 170.54, 40.61, 345.57, 143.47, 242.37, 86.01, 15.04, 186.67, 319.74, 547.8, 1044.63, 1040.28, 43.01, 153.95, 93.7, 62.23, 109.21, 91.74, 3748.94, 350, 672.79, 220.01, 231.88, 153.95, 19523.12, 6372.23, 11741.66, 821.48, 98.7, 443.38, 35.52, 1617.31, 679.46, 727.28, 152.62, 339, 602.09, 269.84, 92.71, 2642.31), ncol = 1, dimnames = list(c("Archaeopteryx", "Jeholornis", "Jinguofortis", "Chongmingia", "Sapeornis", "Confuciusornis_sanctus", "Changchengornis_hengdaoziensis", "Eoconfuciusornis_zhengi", "Confuciusornis_dui", "Yangavis", "Boluochia", "Concornis", "Eoalulavis", "Cathayornis", "Eoenantiornis", "Gobipteryx", "Longipteryx", "Longirostravis", "Neuquenornis", "Pengornis", "Eopengornis", "Chiappeavis", "Parapengornis", "Protopteryx", "Rapaxavis", "Shanweiniao", "Vescornis", "Piscivorenantiornis", "Linyiornis", "Sulcavis", "Bohaiornis", "Longusunguis", "Shenqiornis", "Zhouornis", "Parabohaiornis", "Fortunguavis", "Pterygornis", "Cruralispennia", "Monoenantiornis", "Archaeorhynchus", "Schizooura", "Bellulornis", "Jianchangornis", "Longicrusavis", "Apsaravis", "Hongshanornis", "Archaeornithura", "Parahongshanornis", "Tianyuornis", "Patagopteryx", "Yixianornis", "Piscivoravis", "Iteravis", "Gansus", "Ichthyornis", "Hesperornis", "Baptornis_advenus", "Baptornis_varneri", "Vegavis", "Shangyang", "Mengciusornis", "Mirusavis", "Yanornis", "Similiyanornis", "Abitusavis", "Eogranivora", "Gretcheniao", "Xinghaiornis", "Dingavis", "Qiliania", "Vorona"), "mass"))
mass <- as.data.frame(mass)
mass$mass <- log10(as.numeric(mass$mass))

Next we can generate a PCoA space.

mass.pcoa.data <- Claddis::ordinate_cladistic_matrix(nexus.data)$vectors

Before we go further we should isolate our data to just the taxa we can use (those possessing both all calculable pairwise distances and a mass estimate).

mass.names <- intersect(rownames(mass.pcoa.data), rownames(mass))

And make objects we can use with just these taxa present.

mass <- mass[mass.names, , drop = FALSE]
mass.pcoa.data <- mass.pcoa.data[mass.names, ]
mass.tip.ages <- tip.ages[mass.names, ]
mass.tip.dated.tree <- Claddis::drop_time_tip(tip.dated.tree, setdiff(tip.dated.tree$tip.label, mass.names))
mass.morph.distances <- as.dist(Claddis::calculate_morphological_distances(nexus.data)$distance_matrix[mass.names, mass.names])
mass.mass.distances <- dist(mass[mass.names, "mass"])
attr(mass.mass.distances, "Labels") <- mass.names

Now we can perform various comparisons of mass and morphological variance.

mass.corphy.pco1 <- ape::corphylo(cbind(mass[mass.names, ], mass.pcoa.data[mass.names, 1]), phy = mass.tip.dated.tree)
mass.rlm.pco1 <- MASS::rlm(mass.pcoa.data[mass.names, 1] ~ mass[mass.names, ], as.data.frame(mass.pcoa.data), method = "MM", maxit = 1000)
mass.rlm.pco1 <- lm(mass.pcoa.data[mass.names, 1] ~ mass[mass.names, ], mass.rlm.pco1$model, weights = mass.rlm.pco1$w)
mass.n.pgls <- sapply(1:ncol(mass.pcoa.data), function(idx) {geomorph::procD.pgls(pco ~ mass, mass.tip.dated.tree, data = list(pco = mass.pcoa.data[, idx], mass = mass[mass.names, ]))})
mass.all.pgls <- geomorph::procD.pgls(mass.pcoa.data ~ mass[mass.names, ], mass.tip.dated.tree)
mass.mantel <- ade4::mantel.rtest(mass.morph.distances, mass.mass.distances, nrepet = 9999)
mass.lm.distances <- lm(mass.morph.distances ~ mass.mass.distances)

We can then compile these results in the Brougham and Campione (2020) style.

all_results <- c()
results <- list(pcoa = mass.pcoa.data, corphy_bm_pco1 = mass.corphy.pco1, pco1_rlm = mass.rlm.pco1, pco1_pv = Claddis::ordinate_cladistic_matrix(nexus.data)$values$Rel_corr_eig[1], pco_n_pgls = mass.n.pgls, pco_all_pgls = mass.all.pgls, mantel = mass.mantel, dists_lm = mass.lm.distances, fn = "matrix.nex")
if (!length(all_results)) {
  all_results <- list(results)
} else {
  all_results[length(all_results) + 1] <- list(results)
}
outgroups <- list("matrix.nex" = c("Archaeopteryx"))
all_results <- lapply(all_results, function(x) {x$outgroups <- outgroups[x$fn][[1]]; x})

We can better understand the results through visualisations. Forstly a simpel correlation of body mass and PCo 1.

par(mai = c(0.2, 0.2, 0.2, 0.2), oma = c(4, 4, 0, 0))
stat_table <- data.frame(row.names = c("df", "% PCo1 var.", "Phy. Corr.", "Slope.1", "r2.1", "Slope.2", "r2.2", "Statistic", "r2.3", "pco.all.aov r2"))
idx <- 1
for (result in all_results) {
  plot_bg <- ifelse(rownames(result$pcoa) %in% result$outgroups, "white", "black")
  plot(result$pcoa[, 1] ~ mass[mass.names, ], pch = 21, cex = 0.8, xlab = "", ylab = "", frame.plot = FALSE, col = "black", bg = plot_bg)
  clip(min(mass), max(mass), min(result$pcoa[, 1]), max(result$pcoa[, 1]))
  abline(result$pco1_rlm, col = "red")
  title(chartr("123456789", "abcdefghi", as.character(idx)), adj = 0.05, line = -0.5)
  idx <- idx + 1
  pco1_rlm_stars <- gtools::stars.pval(coef(summary(result$pco1_rlm))[8])[1]
  pco1_r2 <- summary(result$pco1_rlm)$r.squared
  pco1_pgls_stars <- gtools::stars.pval(result$pco_n_pgls[, 1]$aov.table$`Pr(>F)`[1])[1]
  mantel_stars <- gtools::stars.pval(result$mantel$pvalue)[1]
  pco_all_pgls_stars <- gtools::stars.pval(result$pco_all_pgls$aov.table$`Pr(>F)`[1])[1]
  stat_table[result$fn] = c(result$pco1_rlm$df.residual, signif(result$pco1_pv * 100, 2), signif(result$corphy_bm_pco1$cor.matrix[2], 2), paste(signif(result$pco1_rlm$coefficients[2], 2), pco1_rlm_stars), signif(pco1_r2, 2), paste(signif(result$pco_n_pgls[,1]$pgls.coefficients[2], 2), pco1_pgls_stars), signif(result$pco_n_pgls[,1]$aov.table$Rsq[1], 2), paste(signif(result$mantel$obs, 2), mantel_stars), signif(result$mantel$obs**2, 2), paste(signif(result$pco_all_pgls$aov.table$Rsq[1], 2), pco_all_pgls_stars))
}
mtext('Log body mass (kg)', side = 1, outer = TRUE, line = 2)
mtext('PCo1', side = 2, outer = TRUE, line = 2)

We can also compile our results ina stats table.

names(stat_table) <- c("Mesozoic birds")
stat_table <- t(stat_table)
colnames(stat_table)[c(5, 7, 9, 10)] <- "r2"
colnames(stat_table)[c(4, 6)] <- "Slope"
knitr::kable(stat_table, caption = "The results of the robust regression and multivariate phylogenetic GLS for PCo1, and Mantel tests and multivariate phylogenetic GLS of all PCo axes, for discrete character-taxon matrices of Mesozoic birds.", booktabs = TRUE) %>% kableExtra::kable_styling(latex_options = c("striped", "scale_down"), stripe_color = "lightergray") %>% kableExtra::add_header_above(c(" " = 4, "Robust Linear Regression" = 2, "Phylogenetic GLS" = 2, "Mantel test" = 2, "Phylogenetic GLS" = 1))

We can make a barplot of correlation coefficients of body mass and each PCoA axis (accounting for phylogeny).

par(mai = c(0.2, 0.2, 0.2, 0.2), oma = c(4, 4, 0, 0))
idx <- 1
for(result in all_results) {
  stdv <- qnorm(1 - 0.05 / 2) / sqrt(ncol(result$pco_n_pgls))
  coeffs <- apply(result$pco_n_pgls, 2, function (n) {n$aov.table$Rsq[1]})
  barplot(coeffs, col = gray(0.6), ylim = c(0, 1))
  clip(0, ncol(result$pco_n_pgls) * 1.2, min(coeffs), max(coeffs))
  abline(h = stdv, col = "red")
  title(chartr("123456789", "abcdefghi", as.character(idx)), adj = 0.05, line = 0.5)
  idx <- idx + 1
}
mtext('PCoA Axis', side = 1, outer = TRUE, line = 2)
mtext('pGLS correlation coefficient', side = 2, outer = TRUE, line = 2)

We can use our Mantel results to test the signififcance of a mass differences and morphological distances correlation.

par(mai = c(0.2, 0.2, 0.2, 0.2), oma = c(4, 4, 0, 0))
idx <- 1
for (result in all_results) {
  plot(result$mantel$plot$hist, xlim = result$mantel$plot$xlim, main = NULL, col = gray(0.6))
  y0 <- max(result$mantel$plot$hist$counts)
  lines(c(result$mantel$obs, result$mantel$obs), c(y0 / 2, 0))
  points(result$mantel$obs, y0 / 2, pch = 16, cex = 2)
  title(chartr("123456789", "abcdefghi", as.character(idx)), adj = 0.05, line = -0.5)
  idx <- idx + 1
}
mtext('Distance', side = 1, outer = TRUE, line = 2)
mtext('Frequency', side = 2, outer = TRUE, line = 2)

We can also make a plot of morpholoigcal distances versus mass differences.

par(mai = c(0.2, 0.2, 0.2, 0.2), oma = c(4, 4, 0, 0))
idx <- 1
palette <- RColorBrewer::brewer.pal(8, "Set1")
for (result in all_results) {
  dens <- MASS::kde2d(result$dists_lm$model$mass.morph.distances, result$dists_lm$model$mass.mass.distances, n = 100)
  myPal <- colorRampPalette(c("white", palette[c(2, 6, 5, 1)]))
  plot(result$dists_lm$model$mass.morph.distances, result$dists_lm$model$mass.mass.distances, type = "n", frame.plot = FALSE)
  image(dens, col = myPal(256), add = TRUE)
  title(chartr("123456789", "abcdefghi", as.character(idx)), adj = 0.05, line = -0.5)
  clip(min(result$dists_lm$model$mass.morph.distances), max(result$dists_lm$model$mass.morph.distances), min(result$dists_lm$model$mass.mass.distances), max(result$dists_lm$model$mass.mass.distances))
  abline(result$dists_lm)
  idx <- idx + 1
}
mtext('Discrete character distances', side = 1, outer = TRUE, line = 2)
mtext('Body mass distances', side = 2, outer = TRUE, line = 2)

Character exhaustion

We can also calculate character exhaustion values for the full data set and the Enantiornithes and Ornithuromorpha subsets. Note that polymorphisms and uncertainties are not allowed by the function so these are replaced with missing values.

full.exhaustion.results <- paleotree::accioExhaustionCurve(phyloTree = Claddis::drop_time_tip(tip.dated.tree, c("Gallus", "Anas")), charData = apply(nexus.data$matrix_1$matrix, 2, function(x) {if(length(grep("&|/", x)) > 0) {x[grep("&|/", x)] <- NA}; if(any(is.na(x))) x[is.na(x)] <- "?"; if(any(x == "")) x[x == ""] <- "?"; x})[setdiff(rownames(nexus.data$matrix_1$matrix), c("Gallus", "Anas")), ], outgroup = "Dromaeosauridae", missingCharValue = "?", charTypes = gsub("ord", "ordered", nexus.data$matrix_1$ordering))
Enantiornithes.exhaustion.results <- paleotree::accioExhaustionCurve(phyloTree = Claddis::drop_time_tip(tip.dated.tree, setdiff(tip.dated.tree$tip.label, c("Archaeorhynchus", taxon_groups$Enantiornithes))), charData = apply(nexus.data$matrix_1$matrix, 2, function(x) {if(length(grep("&|/", x)) > 0) {x[grep("&|/", x)] <- NA}; if(any(is.na(x))) x[is.na(x)] <- "?"; if(any(x == "")) x[x == ""] <- "?"; x})[c("Archaeorhynchus", taxon_groups$Enantiornithes), ], outgroup = "Archaeorhynchus", missingCharValue = "?", charTypes = gsub("ord", "ordered", nexus.data$matrix_1$ordering))
Ornithuromorpha.exhaustion.results <- paleotree::accioExhaustionCurve(phyloTree = Claddis::drop_time_tip(tip.dated.tree, setdiff(tip.dated.tree$tip.label, c("Protopteryx", taxon_groups$Ornithuromorpha))), charData = apply(nexus.data$matrix_1$matrix, 2, function(x) {if(length(grep("&|/", x)) > 0) {x[grep("&|/", x)] <- NA}; if(any(is.na(x))) x[is.na(x)] <- "?"; if(any(x == "")) x[x == ""] <- "?"; x})[c("Protopteryx", taxon_groups$Ornithuromorpha), ], outgroup = "Protopteryx", missingCharValue = "?", charTypes = gsub("ord", "ordered", nexus.data$matrix_1$ordering))

And visualise them as exhaustion curves.

par(mfrow = c(1, 3))
paleotree::charExhaustPlot(exhaustion_info = full.exhaustion.results, changesType = "charAlt", main = "All taxa")
paleotree::charExhaustPlot(exhaustion_info = Enantiornithes.exhaustion.results, changesType = "charAlt", main = "Enantiornithes")
paleotree::charExhaustPlot(exhaustion_info = Ornithuromorpha.exhaustion.results, changesType = "charAlt", main = "Ornithuromorpha")

Testing for phylogenetic correlation

A potential concern is that our morphological disparity results merely reflect the result of phylogenetic “disparity” between species.

To explore this we can begin by looking at the phylogenetic distances (in millions of years) between our species.

phylo.distances <- ape::cophenetic.phylo(Claddis::drop_time_tip(tip.dated.tree, setdiff(tip.dated.tree$tip.label, rownames(pcoa.data))))
phylo.pcoa.data <- ape::pcoa(phylo.distances, correction = "cailliez")$vectors.cor[rownames(pcoa.data), ]
morph.distances <- Claddis::calculate_morphological_distances(nexus.data)$distance_matrix[rownames(pcoa.data), rownames(pcoa.data)]
morph.pcoa.data <- ape::pcoa(morph.distances, correction = "cailliez")$vectors.cor[rownames(pcoa.data), ]

A simple bivariate plot reveals the issue: our distances seem correlated.

plot(phylo.distances, morph.distances, log = "x", pch = 20, xlab = "Phylogenetic distance (Ma)", ylab = "Morphological distance (arcsine square root of MORD)")

## Warning in xy.coords(x, y, xlabel, ylabel, log): 75 x values <= 0 omitted from logarithmic plot

We also see similar mean paiwrise distances by group…

Claddis::calculate_MPD(list(distance_matrix = morph.distances), taxon_groups)
Claddis::calculate_MPD(list(distance_matrix = phylo.distances), taxon_groups)

…or time bin.

Claddis::calculate_MPD(list(distance_matrix = morph.distances), taxon_groups = Claddis::assign_taxa_to_bins(taxon_ages = pcoa.ages.data, time_bins))
Claddis::calculate_MPD(list(distance_matrix = phylo.distances), taxon_groups = Claddis::assign_taxa_to_bins(taxon_ages = pcoa.ages.data, time_bins))

We can also build dispRity objects as before, for both time bins and groups.

morph.group.disprity <- dispRity::custom.subsets(morph.pcoa.data, group = list(NonOrnithothoraces = match(intersect(NonOrnithothoraces, rownames(morph.pcoa.data)), rownames(morph.pcoa.data)), Enantiornithes = match(intersect(Enantiornithes, rownames(morph.pcoa.data)), rownames(morph.pcoa.data)), Ornithuromorpha = match(intersect(Ornithuromorpha, rownames(morph.pcoa.data)), rownames(morph.pcoa.data))))
phylo.group.disprity <- dispRity::custom.subsets(phylo.pcoa.data, group = list(NonOrnithothoraces = match(intersect(NonOrnithothoraces, rownames(phylo.pcoa.data)), rownames(phylo.pcoa.data)), Enantiornithes = match(intersect(Enantiornithes, rownames(phylo.pcoa.data)), rownames(phylo.pcoa.data)), Ornithuromorpha = match(intersect(Ornithuromorpha, rownames(phylo.pcoa.data)), rownames(phylo.pcoa.data))))
morph.chrono.disprity <- dispRity::chrono.subsets(data = morph.pcoa.data, tree = Claddis::drop_time_tip(tip.dated.tree, setdiff(tip.dated.tree$tip.label, rownames(pcoa.data))), method = "discrete", time = time.bins, inc.nodes = FALSE, FADLAD = disprity.ages.data)
phylo.chrono.disprity <- dispRity::chrono.subsets(data = phylo.pcoa.data, tree = Claddis::drop_time_tip(tip.dated.tree, setdiff(tip.dated.tree$tip.label, rownames(pcoa.data))), method = "discrete", time = time.bins, inc.nodes = FALSE, FADLAD = disprity.ages.data)

And bootstrap them.

morph.group.disprity.boot <- dispRity::boot.matrix(morph.group.disprity, bootstrap = 1000)
phylo.group.disprity.boot <- dispRity::boot.matrix(phylo.group.disprity, bootstrap = 1000)
morph.chrono.disprity.boot <- dispRity::boot.matrix(morph.chrono.disprity, bootstrap = 1000)
phylo.chrono.disprity.boot <- dispRity::boot.matrix(phylo.chrono.disprity, bootstrap = 1000)

And calculate our normal suite of disparity metrics.

morph.group.disprity.boot.sov <- dispRity::dispRity(morph.group.disprity.boot, metric = c(sum, variances))
phylo.group.disprity.boot.sov <- dispRity::dispRity(phylo.group.disprity.boot, metric = c(sum, variances))
morph.chrono.disprity.boot.sov <- dispRity::dispRity(morph.chrono.disprity.boot, metric = c(sum, variances))
phylo.chrono.disprity.boot.sov <- dispRity::dispRity(phylo.chrono.disprity.boot, metric = c(sum, variances))
morph.group.disprity.boot.sor <- dispRity::dispRity(morph.group.disprity.boot, metric = c(sum, ranges))
phylo.group.disprity.boot.sor <- dispRity::dispRity(phylo.group.disprity.boot, metric = c(sum, ranges))
morph.chrono.disprity.boot.sor <- dispRity::dispRity(morph.chrono.disprity.boot, metric = c(sum, ranges))
phylo.chrono.disprity.boot.sor <- dispRity::dispRity(phylo.chrono.disprity.boot, metric = c(sum, ranges))
morph.group.disprity.boot.moc <- dispRity::dispRity(morph.group.disprity.boot, metric = c(median, centroids))
phylo.group.disprity.boot.moc <- dispRity::dispRity(phylo.group.disprity.boot, metric = c(median, centroids))
morph.chrono.disprity.boot.moc <- dispRity::dispRity(morph.chrono.disprity.boot, metric = c(median, centroids))
phylo.chrono.disprity.boot.moc <- dispRity::dispRity(phylo.chrono.disprity.boot, metric = c(median, centroids))

And plot our results. First sum of variance.

par(mfrow = c(2, 2))
plot(morph.group.disprity.boot.sov, main = "Morphology", ylab = "Sum of Variance")
plot(morph.chrono.disprity.boot.sov, main = "Morphology", ylab = "Sum of Variance")
plot(phylo.group.disprity.boot.sov, main = "Phylogeny", ylab = "Sum of Variance")
plot(phylo.chrono.disprity.boot.sov, main = "Phylogeny", ylab = "Sum of Variance")

Sum of ranges.

par(mfrow = c(2, 2))
plot(morph.group.disprity.boot.sor, main = "Morphology", ylab = "Sum of Ranges")
plot(morph.chrono.disprity.boot.sor, main = "Morphology", ylab = "Sum of Ranges")
plot(phylo.group.disprity.boot.sor, main = "Phylogeny", ylab = "Sum of Ranges")
plot(phylo.chrono.disprity.boot.sor, main = "Phylogeny", ylab = "Sum of Ranges")

And centroid medians.

par(mfrow = c(2, 2))
plot(morph.group.disprity.boot.moc, main = "Morphology", ylab = "Centroid Medians")
plot(morph.chrono.disprity.boot.moc, main = "Morphology", ylab = "Centroid Medians")
plot(phylo.group.disprity.boot.moc, main = "Phylogeny", ylab = "Centroid Medians")
plot(phylo.chrono.disprity.boot.moc, main = "Phylogeny", ylab = "Centroid Medians")

rmarkdown::render(input = “~/Documents/Publications/Submitted/Enantiornithes vs Ornithuromorphs - Min/ProjectAncientBeak/markdown/discrete_characters.Rmd”, output_format = “html_document”, output_file = “~/Documents/Publications/Submitted/Enantiornithes vs Ornithuromorphs - Min/ProjectAncientBeak/markdown/discrete_characters.html”)

Bibliography

Brougham, T, and N E Campione. 2020. “Body Size Correlates with Discrete-Character Morphological Proxies.” Paleobiology 46: 304–19.

Brusatte, S L, Benton M J, M Ruta, and G T Lloyd. 2008. “Superiority, Competition, and Opportunism in the Evolutionary Radiation of Dinosaurs.” Science 321: 1485–8.

Foote, M. 1993. “Discordance and Concordance Between Morphological and Taxonomic Diversity.” Paleobiology 19: 185–204.

Gower, J C. 1966. “Some Distance Properties of Latent Root and Vector Methods Used in Multivariate Analysis.” Biometrika 53: 325–38.

———. 1971. “A General Coefficient of Similarity and Some of Its Properties.” Biometrics 27: 857–71.

Hopkins, M J, and K St John. 2018. “A New Family of Dissimilarity Metrics for Discrete Character Matrices That Include Inapplicable Characters and Its Importance for Disparity Studies.” Proceedings of the Royal Society of London B 285: 20181784.

Laurin, M. 2004. “The Evolution of Body Size, Cope’s Rule and the Origin of Amniotes.” Systematic Biology 53: 594–622.

Lloyd, G T. 2016. “Estimating Morphological Diversity and Tempo with Discrete Character-Taxon Matrices: Implementation, Challenges, Progress, and Future Directions.” Biological Journal of the Linnean Society 118: 131–51.

Pyron, R A. 2011. “Divergence Time Estimation Using Fossils as Terminal Taxa and the Origins of Lissamphibia.” Systematic Biology 60: 466–81.

Ronquist, F, S Klopfstein, L Vilhelmsen, S Schulmeister, D L Murray, and A P Pasnitsyn. 2012. “A Total-Evidence Approach to Dating with Fossils, Applied to the Early Radiation of the Hymenoptera.” Systematic Biology 61: 973–99.

Team, R Core. 2019. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org.

Wills, M A, D E G Briggs, and R A Fortey. 1994. “Disparity as an Evolutionary Index: A Comparison of Cambrian and Recent Arthropods.” Paleobiology 20: 93–130.